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ABSTRACT 

It is well-recognized that the presence of magnetic fields will lead to 
anisotropic energy cascade and dissipation of astrophysical turbulence. With 
the diffusion approximation and linear dissipation rates, we study the cascade 
and damping of Alfven-cyclotron fluctuations in solar plasmas numerically for 
two diagonal diffusion tensors, one (isotropic) with identical components for the 
parallel and perpendicular directions (with respect to the magnetic field) and 
one with different components (non-isotropic). It is found that for the isotropic 
case the steady-state turbulence spectra are nearly isotropic in the inertial range 
and can be fitted by a single power- law function with a spectral index of —3/2, 
similar to the Iroshnikov-Kraichnan phenomenology, while for the non-isotropic 
case the spectra vary greatly with the direction of propagation. The energy fluxes 
in both cases are much higher in the perpendicular direction than in the parallel 
direction due to the angular dependence (or inhomogeneity) of the components. 
In addition, beyond the MHD regime the kinetic effects make the spectrum softer 
at higher wavenumbers. In the dissipation range the turbulence spectrum cuts off 
at the wavenumber, where the damping rate becomes comparable to the cascade 
rate, and the cutoff wavenumber changes with the wave propagation direction. 
The angle averaged turbulence spectrum of the isotropic model resembles a bro- 
ken power-law, which cuts off at the maximum of the cutoff wavenumbers or 
the 4 He cyclotron frequency. Taking into account the Doppler effects, the model 
naturally reproduces the broken power-law turbulence spectra observed in the 
solar wind and predicts that a higher break frequency always comes along with 
a softer dissipation range spectrum that may be caused by the increase of the 
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turbulence intensity, the reciprocal of the plasma /3 P , and/or the angle between 
the solar wind velocity and the mean magnetic field. These predictions can be 
tested by detailed comparisons with more accurate observations. 

Subject headings: MHD — plasmas — solar wind — turbulence — waves 



Introduction 



Turbulence is ubiquitous in t he universe and pla y s important roles in our understand - 
ing of many natural phenomena (IKolmogorov I Il94ll ; llroshnikov I Il963l : iKraichnan I Il965l ) . 
It occurs in highly non-equilibrium systems, where the microscopic viscous and/or resistive 
dissipation processes cannot effectively convert the free energy into the internal energy of the 
fluid. Such systems usually have very high Reynolds numbers and/or cover huge spatial and 
temporal scales, so that the free energy is stored in the large scale motions and/or magnetic 
fields. For plasmas, the free energy initially may also be stored in non-equilibrium distribu- 
tions of charged particles. Turbulence is generated through a variety of instabilities related 
either to the large scale magneto-hydrodynamic (MHD ) processes or the microsco pic col- 
lective plasma effects and/or plasma physics processes (IBorovsky fc Funsten 1 120031 ) . These 
aspects have been extensively investigated with the MHD and/or plasma physics theories. 
In astrophysics, most turbulence is carried by magnetized plasmas. It is responsible for dis- 
tributing energies among different components of the plasmas, which may result in distinct 
emission characteristics or other observable features. Observations of these radiations can 



be used to study the corresponding astrophysical sources [e.g., iLiu et al. I (12004 . 120061 ) on 
3 He rich impulsive solar energetic particle events]. 

Plasma waves or turbulence once generated (usually on large scales comparable to the 
dynamical scales given the large amount of energy inferred from astrophysical observations) 
undergo two types of interactions: wave-wave and wave-particle. The former dominates at 
large scales, where the Reynolds (ordinary or magnetic or both) number is large, and re- 
sults in an almost dissipationless cascade primarily to smaller scales and higher values of the 
wavenumber k. The interactions of waves with the more numerous low energy particles grad- 
ually become more important as the cascade reaches short enough scales (and the Reynolds 
number approaches unity and/or the kinetic effects set in) and damp the waves. This results 
in plasma heating and/or particle acceleration. The wave-particle interactions also deter- 
mine the spatial diffusion and energizing of high energy par ticles in collisionless turbulent 



plasmas and are essent ial elements of the cosmic ray theory (IGinzburg fe Syrovatskii 111969 



Yan &: Lazarian Il2002l ). The spectral range from the turbulence generation scale ko to the 
damping scale /c max is called the inertial range, where the angle averaged energy flux is inde- 
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pendent of k. Beyond this range, there are strong couplings between the charged background 
particles and turbulent motions. For collisional plasmas, where the Coulomb collision mean 
free path is much shorter than other relevant spatial scales, the energy from the turbulence 
dissipation is thermalized quickly and the result is plasma heating. Much astrophysical tur- 
bulence, however, is carried by collisionless plasmas, where the Coulomb mean free path is 
long. Although it is generally accepted that the particle distributions in these plasmas are 
determined by couplings of the (charged) particles with the turbulent electro-magnetic fluc- 
tuations, the details of these interactions are not well understood. The wave damping, which 
usually occurs on the smallest spatial scales, has been an essent ial part of plasma physics 



theori e s for collisionless bu t presu mably thermal plasmas [see e.g. , iBraginskii I (119651 ) ; lAndre 



( 119851 ); iGary &: Borovsky I (12004)] . There i s littl e work on damping by nonthermal particle 



distributions. Recently, iPetrosian et al. I (120061 ) derived the transit-time damping (TTD) 
rate of fast-mode waves by an arbitrary energy distribution of electrons and protons under 
typical solar-flare conditions. Large scale waves in this case can be damped by high en- 
ergy particles through resonant wave-particle couplings. The studies of the turbulent energy 
dissipation in general and the particle acceleration from a low energy background to high 
energies in particular require a more complete (beyond the MHD regime) treatment of the 
wave cascade and damping. 

The cascade process has been an essential element in all kinds of turbulence studies. The 
highly nonlinear nature of turbulence makes this a very challenging aspect of the problem. 
Nevertheless, the energy transfer of isotropic incompressible hydrodynamical turbulence, 
the simplest form of turbulence, has been described reasonably well with the Kolmogorov 
phenomenology that assumes a scale-independent self-similar cascade process, leading to 
the famous Kolmogorov power-law spectrum with a ( one dimensional isotropic) spectral 
index of —5/3 in the inertial range (IKolmogorov Ill94ll ). Most astrophysical plasmas carry 
magnetic fields so that large scale MHD waves and small scale kinetic plasma fluctuations 
may be excit ed. The wave prop agat ion effects on the cascade of MHD turbulence were first 
discussed by llroshnikov I (119631 ) and iKraichnan I (119651 ). Although MHD effects are expected 
to introduce anisotropy to the system that can affect the couplings between the turbulence 
and background particles significantly, these preliminary investigations assumed that the 
turbulence is isotropic a nd reduced the cascade to a one-dimensional (ID) problem, yielding 
a spectral index of —3/2 (IZhou fc Matthaeus Ill990l ). This is called the Iroshnikov-Kraichnan 
(IK) phenomenology. Since both the Kolmogorov and IK cascade timescales (oc k~ 2 ^ 3 and 
fc _1//2 , respectively) decline with increasing k (or decreasing spatial scales) more slowly than 
the periods of MHD waves (oc A; -1 ), the turbulence may be better described as spectra 
of waves at higher values of k's (smaller spatial scales), which play a critical role in the 
energizing of low energy background particles. However, at such high values of fc's one may 
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be stepping beyond the MHD regime and must use more complex dispersion r elations to take 



into account the kinetic effects and the anisotropy of the turbulence properly (ILeamon et al 



20001 ; iMarkovskii et al. 112000 ). 



To study the anisotropic effects induced by the presence of large scale mag netic fields, 



MHD turbulence has been treated as an ensemble of linear wave modes [e.g., ICho et al. 



( 120031 )]. Its nonlinear nature is revealed in the wave- wave couplings. Although this quasi- 
linear treatment may not be valid for strong turbulence, it is certainly a good approximation 
when the magnetic field fluctuations are much smaller than th e large scale fi eld and the wave 
periods are much shorter than the eddy turnover timescales (IGaltier 1120061 ). Significant in- 
sight of propert ies of Alfven and magnetosonic turbulence has been obta i ned recently through 



this approach (jSridhar &: Goldreich I Il994l ; iGoldreich fc Sridhar I Il995l ; IGaltier et al. 1 12000 



Chandran 20051 ). However, the component of Alfven turbulence exci tations nearly perpen- 



dicu 



ar to a large scale magnetic field is dominated by nonlinear effects ^Montgomery fc Turner 



1981 ). This ind icates inherent limitations of treating MHD t urbul ence as a spectrum of waves 



Galtier II2006I ). As shown by IMontgomery fc Matthaeus I (119951 ). linear MHD wave modes 
do not give a complete description of turbulence excitations and their couplings, and treat- 
ing MH D turbulence as an ensemb le of linear wave modes may miss some critical nonlinear 
effects (jNg fc Bhattacharjee 1997 ). This is more evident in hydrodynamic simulations of 
decaying turbulence, wher e the compressible wave co mponent accounts for about 10% of the 



turbulence kinetic energy (jPorter et al. 



( IMatthaeus et al. Ill996l ; IVestuto et al. 



19981 . 1999 ). MHD simulations give similar results 
20031 ) . These studies suggest that treating turbu- 



lence as a spectrum of linear wave modes may not describe the dominant energy component 
properly, which also reflects the highly non-linear nature of a fully developed turbulence, es- 
pecially o n large scales, where the way e periods may be longer than the turbulence cascade 
timescale (IBorovsky fc Funsten 1120031 ). 



Based on the dominance of nonlinear or wave-propagation effects, lOughton et al. I (120061 ) 
recently separated the Alfven turbulence into two interacting parts: quasi-two-dimensional 
(2D) and wave-like fluctuations. The quasi-2D component characterizes the nearly per- 
pendicular excitations and may be described with a quasi-2D Kolmogorov phenomenology. 
The wave-like fluctuations may be described with an IK phenomenology with propagation 
direction dependent wave periods. The cascade of turbulence then also depends on the 
assumed couplings of these two components. It is interesting to note that the two com- 
ponents are separated by the critical balance between linear wave periods and nonlinear 
eddy turnover timescales, a nd the incompressible strong Alfven turbulence described by 
Goldreich fc Sridhar I (119951 ) appears to be appropriate for the quasi-2D component that 
is dominated by nonlinear effects. Due to the suppression of cascade by wave propaga- 
tion effects, Alfven turbulence with wave periods proportional to the parallel component of 
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wavevectors cascades preferentially in the direction perpendicular to the mean magnetic field. 
The wave-like fluctuations need longer time to develop, which may explain the exponential 
cutoff of the power spectrum i n the direction pa rallel to the mean large scale magnetic field 



observed in MHD simulations (jCho et al. Il2002l ). 



Advances in computational power and numerical algorithms over the past few decades 
have made numerical si mulations one of the important tools for quantitative investigations of 



magnetized turbulence (jShebalin et al. Hl983l ; lMatthaeus et al. lll998l ; ICho et al. 1120021 . 12003 



Cho fc Lazarian II2003I ). However, turbulence usually covers a huge dynamical range from the 
macroscopic scales of turbulence generation to the microscopic dissipation scales. Current 
simulations have a dynamical range of a few hundreds to a few thousands and have not been 
able to give a complete description of energy flows in magnetized turbulence, limiting their 
astrophysical applications. Moreover, most of these studies are limited to the MHD regime, 
where the background particles are strongly coupled with each other and can be treated 
together as a single fluid. The electron magnetohydrodynamics (EMHD) treats electrons 
and ions as two fluids. It is only applicable in a narrow frequency range bet ween the electron 



and i on gyro-frequencies, w here the whistler dispersion relation is valid (iBiskamp et al 



19991 ; iPetrosian fc Liu 1 12004 ) . Recently Galtier (2006) showed that the transition from the 
MHD to the EMHD may be studied with the incompressible Hall MHD for weak turbulence. 
However, a quantitative weak turbulence theory is still to be developed, and it is not obvious 
how a complete theory may cover both the weak and strong turbulence regimes. To address 
the heating of background particles by turbulence, one usually assumes that most of the 
background particles reach thermal distributions and arbitrarily extrapolates the turbulence 
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( 120061 )]. The particle-in-cell and kineti c simulations can usually study processes n ear the 



gyro-scales of the background particles (IGary. Saito. fc Li 1120081 ; iHowes et al. II2008I ). 



Observations of solar corona, solar flares, solar winds, and space plasmas, on the other 
hand, demand a detailed study of the turbulence evolution over a large dynamical range. 
There are no direct observations of turbulence in solar corona and solar flares, but the ob- 
served high Reynolds numbers and models for particle acceleration demand existence of tur- 



bulent plasma waves (IBrown fc Melrose 1119771 ; iMiller fc Roberts 1 1 1 9 9 5l : iFletcher fc Hudson 



20081 : iKrucker et al. Il2008bh. In situ observations of space plasmas reveal evidence of wave 



particle interactions (IParks et al. 1120061 ; IChen et al. 1120071 ). turbulence generation and dissi- 
pation, and the d ominance of the turbulence (in the Earth's plasma sheet) by eddies instead 
of plasma waves (IBorovsky fc Funsten 1 120031 ) . In solar winds, direct measurements of the 
turbulence spectrum over more than 6 orders of magnitude in frequency challenge all ex- 



isting theoretical models (IDenskat et al. Ill983l ). Indeed, strong anisotropics are observed 
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in almost all the relevant spatial scales (Hamilton et al. 1 120081 ; lOsman Horbury 1 12007 



Dasso et al. 1120051 ; iMatthaeus et al. Ill990l ). While the turbulence anisotropy at large scales, 



which is likely related to t he turbulence gene ration mechanism, appears to be correlated 



with the solar wind speed (IDasso et al. II2005I ). at smaller scales the quasi-2D component 



always dominates, whi ch should be attributed to the dynamics of the turbulence cascade 



( jHamilton et al. II2008I ). At even higher frequencies, there is a spectral break, which has 



(ILi et al. 112001 



Stawicki et al. 


200ll: 


Galtier 


2006: 


Howes et al. 
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tion range and kinetic effect dominant regime, background particles with different charge to 
mass ratios interact w i th the electromagnetic fluct uations quite differently (jPetrosian fc Liu ' 



2004 ; iLiu et al. 1120041 . 120061 ; iLeamon et al. Ill998l ). This results in complicated wave disper- 



sion relations, namely the dependence of the wave frequency on the wavevector ((Andre 
19851 ). and certainly affects the turbulence cascade. The couplings among wave s and par- 



ticles are also complicated (IStix Ill962l ; IXie fc Ofman 112004 ; ISaito fc Gary 1 120071 ). Thus for 
practical applications one requires a more readily solvable approach and algorithm, which can 
utilize the basic physics learned from detailed simulations in an approximate but numerically 
effective manner. 

The diffusion approximation for the power spectrum in the wavevector space has been 
a very efficient tool to study the turb ulence cascade and dissipation over a large dynamical 
range (IZakharov fc Kuznetsov Ill978l ). Although the diffusion equatio n may not be derived 
from the fluid dynamical equations, especially for weak turbulence ( INg fc Bhattacharjee 



19971 ; iGaltier et al. 1 120001 ) . it has several achievements in quantitative studies of the MHD 
turbulence. The ID models not only address the transition o f turbulence from the large scale 
Kolmogorov phenomenology to the small scale IK cascade (IZhou fc Matthaeus Ill990l ). but 
also are use d to study the accelerat i on and heating of background particles by magnetized 
turbule nce (IMiller fc Roberts Ill995l ; iMiller et al. 1119961 ) and the damping of waves at small 
scales (ILi et al. I l200ll ). These studies have deepened our understanding of energy release 
processes during solar flares significantly. However, they encounter difficulties in reproducing 



the observed broken p ower-law spectrum of solar wind magnetic fluctuations (ILi et al. 112001 
Stawicki et al. Il200ll ). Due to the nonlinear nature of these ID diffusion models, the tur- 
bulence spectrum cuts off sharply at the wavenumber, where the cascade timescale becomes 
comparable to the wave damping time. 

Studies of weak Alfven turbulence with the quasi -linear treatm ent of electromagnetic 
fluctuations have shown that the cascade is anisotropic (IGaltier 1120061 ) . It is also well-known 
that damping rates of different plasma wave modes by thermal ba c kgrou n d particles are very 



sensit iv e to the wave propagat ion direction [see e.g., iGinzburg I (Il96ll ); iGary Borovsky 



( 120041 ); iPetrosian et al. I (120061 )]. These anisotropies are critical to study the energy dis- 
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sipation through magnetized turbulence in a collisionless thermal plasma. Thus, 2D dif- 
fusion models are needed to advance our understandi ng of magnetized turbulence and ad - 
dress the difficulties encountered with the ID models. ICranmer fc van Ballegooijen I (120031 ) 
have shown recently that the heating of the background particles by Alfven turbulence are 
very sensitive to the 2D an d kinetic effects. To recover the critical balance proposed by 
Goldreich fc Sridhar I (119951 ) and partially take into account the kinetic effects, they con- 
structed a complicated advection-diffusion equation for the power spectrum with three di 



mensionless c oefficients. A similar quasi-2D model was proposed recently by iHowes et al. 



( 120071 . l2008al ) to explain the broken power-law character of the power spectrum of the tur- 
bulent magnetic field fluctuations in the solar wind. 

In this paper, we study the general characteristics of the nonlinear 2D diffusion model 
with the kinetic and damping effects of Alfven-cyclotron fluctuations treated self-consistently. 
In principle, all plasma mode branches, such as Alfven, fast, and slow modes, and their 
extension beyond the MHD regime at small scales, should be included to have a complete 
description of the cascade and damping of magnetized turbulence. One then ne eds to solve 
a set of coupled diffusion equations for each wave mode branch (lAndre I Il985l ). This is a 
quite challenging task because the results will depend on b oth the interactions wit hin each 
wave mode branch and couplings among different branches. ICho fc Lazarian I (20031) showed 
that couplings among different mode branches are usually weak, and Luo fc Melrose ( bo06 ) 
showed that these couplings decrease toward small spatial scales. Ichandran" ( 2005 ). on 
the other hand, found that the couplings between Alfven and fast modes are significant in 
the direction parallel to the large scale magnetic field, where the frequencies of the Alfven 
and fast modes are comparable. The couplings among different branches therefore can be 
separated from other processes near the dissipation range, where the kinetic and damping 
effects dominate. 

We use the exact dispersion relation for a cold plasma (to go beyond the MHD regime), 
which is a good approximation for the more general dispersion relation of a collisionless ther- 
mal plasma. The diffusion coefficients are constructed as a function of the wave frequency, 
phase or group velocities so that the kinetic effects are treated self-consistently. In general, 
the linear Vlasov equation is used t o derive the therm al damping rates of the waves. Here 
we use the so-called WHAMP code ( Ronnmark 1982 ) for calculation of the damping rateE 
The diffusion model a ctually does not distinguish the quasi-2D and wave-like fluctuations 
(jOughton et al. 1120061 ). The balance between the eddy turnover and wave propagation is 



1 Since the damping rate increases sharply with the increase of the wavenumber, the quasi- linear treatment 
is expected to give a good approximation of the damping even for the nonlinear effect dominated quasi-2D 
component. 
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revealed in the diffusion tensor. 

We apply our formalism (for cascade and damping) to the Alfven- cyclotron branch, 
which has been studied extensively due to its simplicity and its prevalence in magnetized 
turbulence. It is chosen here to facilitate better comparisons of our models with previous 
studies. In § El we discuss how the turbulence cascade can be studied using the diffusion 
approximation, which reduces the turbulence evolution to a nonlinear 2D diffusion problem. 
We also investigate different forms for the diffusion tensor and include wave damping rate. 
For the simplest cases, where the diffusion tensor is only a function of time, wavevector, 
and local turbulence power in the wavevector space, the nonlinear diffusion equation can be 
solved numerically to obtain the power spectrum from the MHD region to the 4 He gyro- 
frequency, where the Alfven dispersion surface cuts off. To better understand the cascade 
models, we first present results for the Alfven turbulence § [21 where the dispersion relation 
for Alfven waves is used. The exact dispersion relation and the cascade beyond the MHD 
regime are discussed in § [U The damping effects are investigated in § [51 In § [6] we apply 
models to solar wind magnetic fluctuations and show that a locally isotropic diffusion model 
can naturally explain the observations. We discuss the implication of this theory on the study 
of plasma heating and particle acceleration by magnetized turbulence and future work, and 
draw conclusions in § [71 



2. Kinetic Equation 



Following the approach of IZhou fc Matthaeus I (I1990I ). we assume local interactions in 
the wavevector k space so that the evolution of the turbulence power spectrum W(k, t) only 
depends on its properties near k at the time t: 



dW(k,t) 
dt 



Q(M) + 



_d_ 



Dy(k,t)A W (k,t) 



-r(k,t)w(k,t) 



^~csc(k, t) 



Here i,j indicate the three bases of the wave vector k, and the terms on the right-hand- 
side represent the turbulence generation, cascade (through diffusion), damping, and leakage 
processes, respectively. The rates of turbulence generation (Q), damping (r), and leakage 
(T~c) are mostly determined by the properties of the background plasma, whose variation 
timescales, comparable to the large scale dynamical time, are usually much longer than the 
turbulence evolution timescales so that the time- dependence of these rates can be ignored. 
The diffusion tensor -Djj(k, t) can be constructed with the nonlinear timescale TNL(k, i) as- 
sociated with the eddy turnover, the wave crossing time rw(k), and the triple correlation 



time T3 - 
A, ~ k 2 



(' 



NL 



w 



The turbulence cascade rate then is defined as r c J. ~ r 



NL T 3- 



and 



T„, 
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More generally, one may also include an advection term to the above equation so that 
the energy flux in the wavevector space 



/i(M) 



-Dy(M)A>V(M) 



U, 



(M)W(M) 



(2) 



where it, is a velocity field in the wavevector space (ICranmer &: van Ballegooijen II2003I ). This 
velocity field can be decomposed into a solenoidal u s and a compressible component u c . The 
former is divergence-free and can be attributed to the anti-symmetric part of the diffusion 
tensor £^-(k, t). The latter is driven by an external potential in the k space $(k, t), i.e., 
u C i = d^/dki. The highly non-linear phenomenon of intermittence can also affect the energy 
diffusion of turbulence. However, int ermittence likely results from the non-local nature of 
interactions in the wavevector space (jBoldyrev 1 120021 ; ICho et al. Il2002l ). One may partially 
take into account the e ffects of this process on the turb ulence energy transfer with a nonlinear 
damping rate [see e.g., iBeresnyak fc Lazarian I (120081 )]. In general, the diffusion tensor is a 
functional of W. For local interactions in the k space, it is reduced to a function of k. In what 
follows, we will focus on the relatively simple case with u c = 0, and Dij symmetric and only 
functions of k and t, and ignore the cascade and dissipation associated with intermittence. 

Without excitation of waves, i.e., r w = oo, the turbulence cascade is well described with 
the Kolmogorov phenomenology, which assumes an isotropic incompressible fluid with local 
interactions of eddies in the wavenumber space so that tnl = (WA; 5 ) -1 / 2 and = Syk 2 /tnl, 
where Sy is the Kronecker's delta function. For strong turbulence with r NL (k) <C T W (k), 
nonlinear processes dominate and one retains the Kolmogorov phenomenology. In the weak 
turbulence limit, tnl (k) r\y(k) and the turbulence may be described as a spectrum of 
linear wave modes. For fast mode waves in a cold or low beta plasma (i.e., when the Alfven 
speed va Cs, the sound speed), the wave frequency u;(k) ~ v&k, and rw = cu(k) -1 <C 
tnl = (WA; 5 ) -1 / 2 , = 5ijkr^l/vA, and one recovers the IK phenomenology^ For Alfven 
waves, u;(k) = i>a&||, where k\\ is the component of the wavevector parallel to the large scale 
magnetic field, wave packets propaga te along the mean m agnetic field. Since only oppositely 



directed wave packets can interact (IGaltier et al. 1 120001 ). rw oc k,, . Although the weak 



turbulence approximation may be valid for parallel propagating waves, nonlinear processes 
always dominate in the perpendicular direction. Therefore the quasi-linear theory for linear 
wave modes does not give a complete description of the Alfvenic turbulence fluctuations 
ffoughton et al. Il2006h . 



We are interested in extending the turbulence cascade study beyond the MHD regime 
by using the general dispersion relation for the Alfven-cyclotron fluctuations and compress- 



2 Note that for sound waves in a fluid, w(k) = Csk, one obtains results similar to fast mode waves. 



-10- 



ible plasma waves. The quasi-linear approaches for MHD wave interactions then become 
extremely complex if not completely unsolvable. On the other hand, the diffusion approxi- 
mation can provide simple solutions with a prescribed diffusion tensor and may help us to 
uncover processes beyond the MHD regime. In the following, we will ignore the leakage term 
in equation (TTJ and treat the source term as a ^-function at certain large scale. Thus, for the 
investigation of the evolution of the turbulence power spectrum we only need proper choices 
of the diffusion tensor Dij and the damping rate T. 



2.1. Diffusion Tensor 

We adopt the recipe given by Zhou & Matthaeus (1990) for the ID cascade of Alfven 
waves caused by the wave-wave resonances and eddy turnovers and generalize it to waves 
propagating in arbitrary directions with a propagating direction dependent diffusion tensor. 
The eddy turnover time can be estimated with tnl = (vkk)^ 1 , where the eddy velocity § 

v k ~ (Wk 3 ) 1/2 so that r NL = (WP)- 1/2 . (3) 

For the wave crossing time, one must consider the anisotropy of the Alfven-cyclotron disper- 
sion. Wave packets with a size of ~ A;^ 1 cross each other at the Alfven speed. We then have 
r w ~ (fA^ii) -1 , and the wave propagation effect reduces the triple correlation time from tnl 
to r 3 and gives a cascade rate r c ~^ = t^ L J {j N \ + r^ 1 )!^ Then the simplest (locally isotropic) 
diffusion tensor may be constructed as: 

C_ 2 _ x C Wk 7 C Wk 6 

ij - i3 4n Tcas i3 47r (Wk?y/ 2 k + v A kcos9 'Un v k + v A cos9 1 ] 

where C is a dimensionless scaling constant corresponding to the Kolmogorov constant 
for hydrodynamic turbulence, and 9 is the angle between k and the mean magnetic field: 
k\\ = k cos 9. In what follows we set C = l|3 This diffusion tensor is locally isotropic in 
the wavevector space implying that the energy of a wave packet at k can spread into its 



3 Here we have ignored the potential anisotropy of the turbulence. In a more accurate treatment, one also 
needs to distinguish contributions to the turbulence power spectrum fro m the velocity and magnetic fields. 



The eddy velocity should only depend on the turbulence kinetic energy ([Galtier 1120061 ). 

4 Note that, following Kraichnan's argument (1965) for an isotropic Alfven wave turbulence, Zhou & 
Matthaeus (1990) obtained rw = (t>A&) -1 for their ID diffusion model, which is appropriate for the isotropic 
sound and fast mode waves. 



5 C = (3/ll)(47r/C*o) 3/2 ~ 5.89, where C 
the eddy speed v k = (47rWfc 3 ) 1 / 2 , C ~ 1.66. 



1.62 is the Kolmogorov constant (jYeung fc Zhou Ill997f ). If 
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neighborhood in the k space without any preferred directions. The wave propagation effect 
just reduces the diffusion coefficient. This reduction is higher for waves with higher phase 
speeds. 

In more general cases, the diffusion tensor may not be locally isotropic in the wave 
vector space due to resonance interactionsfl One may construct the parallel and per- 
pendicular components of the diffusion tensor differently. In the perpendicular direction, 
the turbulence cascade is determined by the eddy turnover time only. We have D± t ± = 
^ 2r cas i = ^ 2t nl- The wave propagation effect reduces the cascade rate in the parallel direc- 
tion r^ii = Tnl/( t nl + r w )> an d t ne corresponding diffusion coefficient can be written as 
Du II = k 2 r~ ,,. Thus, we obtain an anisotropic diffusion tensor 



1 k 2 



tnl/Vw) 1 
1 



(5) 



13 4vr r NL 

There are also two choices for the wave interaction timescale. We can choose 

r w x = k.v gr , or w(k) , (6) 

where v gr is the wave group velocity. In the MHD regime they both give rw = kv a cos 9 
used in equation (J3J). However, as mentioned above we are often interested in acceleration 
of low energy particles, which requires cascade to high values of k regime. In this case we 
should use the exact dispersion relation, going beyond the MHD regime (k > Q p /va), where 
the above timescales are different. In what follows we will present results for both tensors. 
A brief review of the exact dispersion relation and its decomposition into different modes is 
presented in Appendix A. Figure d3] gives the u — k relation for different modes at various 
angles of propagation with respect to the large scale magnetic field. 

So far we have only used the dispersion relation of waves to construct the diffusion 
tensor. Physically the polarization state (including the compressibility) of the waves may 
also affect the wave diffusion. An equation for the magnetic helicity spectrum is also required 



(IGaltier II2006I ). These may become more significant when we consider the couplings between 



different wave branches. They may also explain the weak couplings between incompressibl e 



Alfven wave modes and the compressible fast and slow mode waves ( jCho &: Lazarian 1120031 ) 



and the relatively strong couplings between parallel propagating Alfven and fast mode waves 



(IChandran II2005I ). Instead of exploring more possibilities on the forms of the diffusion tensor, 
we defer to observations of turbulence spectrum in solar winds to point us the right direction. 
The detailed discussions on this is given in § [6] below. 



3 Non-resonance interactions presumably lead to an isotropy. 
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2.2. Damping Rate 



MHD turbulence in collisionless plasmas is subject to damping caused by both res- 
onance and non- resonance processes. The latter may be associated wit h current sheets 



and is not well understood (jLeamon et al. I l2000l; iMarkovskii et a 



likely dominates the damping of plasma waves (jCho &: Lazarian 



2006), and the former 
20061 ). It is well-known 



that charged particles can extract energy from plasma waves through cyclotron resonances, 
transit-time damping (TTD), and Landau damping (LD) processes. Cyclotron resonances 
by low energy background particles can only lead to damping of a high frequency wave with 
the wave frequency comparable to the particle gyro-frequency. Magnetic field fluctuations 
in the direction of the mean magnetic field are required for the TTD processes to operate. 
The TTD conserves the magnetic moment of the particles. Efficient pitch a ngle scatterin g 
of the particles is also needed for this process to energize particles effectively (IMiller 1119971 ). 
The LD is related to electric fiel d fluctuations along the mean magnetic field and therefore 
can be very efficient (jStix Ill962l ). Alfven waves in the MHD regime don't have electric and 
magnetic field components along the mean magnetic field and therefore are only subject to 
cyclotron damping by relatively high energy particles, whose gyro-frequency is relatively low. 
At small scales and/or high frequencies, where the kinetic effects become important, there 
are electric and magnetic field fluctuations along the mean magnetic field, and the TTD and 
LD processes can damp the Alfven wave branch. 

The damping of waves energizes the plasma either by heating it or by production of 
nonthermal tails. The properties of the turbulence and background plasma determine the 
energy partition into these two channels. Therefore, in general, both the hot plasma and 
non-thermal particles contribute to the damping rate. Unfortunately there is very little work 
on determin ation of the damp i ng rat e by nont hermal particles. T he only result we are aware 
of is that by iPetrosian et al. I (120061 ) [see also lYan et al. I (120081 )] on the TTD of fast mode 
waves in low beta plasmas; (3 P < 1. This was done by using the fact that the rate of particle 



acceleration and the dam ping rate are related via the energy conservation (IMiller fc Roberts 



19951 ; IMiller et al. Ill996l ). It was shown that for typical solar flare conditions, except for very 
large flares, the damping is usually dominated by thermal particles. We believe that the 
same is also true in the solar wind, where we apply our formalism in this paper, as indicated 
by the relatively low energy content of non-thermal particles. Thus, in what follows, we 
include damping by thermal particles only, in which case the damping rate as a function of 
wavevector k is given by the imaginary part of the frequency obtained from the dispersion 
relation for a warm (T ^ 0) plasma. 

For thermal background particles with isotropic pitch angle distributions, depending on 
the background magnetic field, particle abundances and temperatures, the different damping 
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processes o f the Alfven wave branch mentioned above dominate in different wave propagating 
directions (IGary fc Borovsky I I2004J ) . At low values of k\\ < Q a / V A, where Q a is the non- 
relativistic gyro-frequency of the a-particle, these waves are called the kinetic Alfven waves 
(KAW) near k±p p ~ 1, where p p is the Larmor radius of the thermal background protons . 



They play important roles in the heating of background particles ( jLeamon et al. 1 11999 



Cranmer fc van Ballegooijen II2005I ). The quasi-parallel propagating Alfven waves have very 
weak electric and magnetic field fluctuations along the large scale magnetic field and therefore 
are not subject to the TTD and LD. These waves play important roles in the selective 
acceleration of backg round ions with different gyro-frequencies to high energies th rough 



cyclotron resonances (IReames et al. 1 1 1994 iMason et al. 



20021 : iLiu et al. 112004 12009 ) . 



Swanson I (119891 ) gives a simplified expression for the damping rate by electron cyclotron 
resonances of parallel propagating waves under the weak damping approximation. This result 
can be easily generalized to calculate the cyclotron damping rate by all background particle 
species: 

Ea/vto;^ r fu — e.,fL 
. — exp 



0J r 



0J r kv 



t.s 



kvt 



2k 2 

UJ 2 ^ 



2uo — e,=fL 



Here subscripts i and r refer to the imaginary and real parts, s stands for different particle 
species and u p , fl, v t , and e are the corresponding particle plasma frequency, non-relativistic 
gyro-frequency, thermal speed, and sign of the charge, respectively. For a given power 
spectrum of quasi-parallel propagating waves, this equation can be used to study the relative 
acceleration of different background particle species. 

We use the Waves in Homogeneous Anisotropic Multitemperature Plasma (WHAMP) 
code lRonnmark I (Il982l ) to obtain the thermal damping rate numerically. When T = Ui ~ u r , 
the weak damping approximation is invalid and the code fails to converge. We use a power- 
law extrapolation with respect to the wavevector to approximate the damping rate in the 
strong dissipation regime. This does not affect our result because the turbulence spectra 
usually cut off at the wavevectors where the damping rate is still much lower than the wave 
frequency. For further details see Appendix B and Figures [5] and [6] below. 



3. Cascade of Alfven Turbulence 



In this section we consider the cascade in the inertial range ignoring damping and 
leakage. To compare with previous studies, we first consider Alfven turbulence in the MHD 
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regime with r^ 1 = u;(k) = k\\vx- In the strong turbulence limit one has v\ va and 
tnl <C r-w- As a result the diffusion tensors of equations (jlj and (J5]) lead to the isotropic 
Kolmogorov cascade. However, in the weak turbulence limit, i.e., <C v^, they are different 
from the isotropic IK phenomenology due to the dependence of the Alfven wave crossing 
time on the angle of propagation 9. Since the Alfven wave crossing time t\v is much longer 
in the perpendicular direction than in the parallel direction, the wave cascade rate is high in 
the perpendicular direction and decreases dramatically with the increase of h\. We solve the 
2D time-dependent wave kinetic equation (CQ) with the alternative direction implicit scheme 
on a uniform log-log grid and with 

Q w (k,t) = F <f(k-ko)e(t), 

k u = k ±j0 = k = 1.6 x lO~ 8 tt p /v A , F = 9.55 x KT 10 ^^, (7) 

where G is the Heaviside step function and S is the 3D Dirac-5-function. A reflective bound- 
ary condition is used at the large scale boundary feii o, &_l,o- I n this section, we ignore the 
damping and leakage and allow the turbulence to escape freely at the high end boundary 
of &;[] Because no damping is included, the time unit is set by the proton gyro-frequency 
Q p . The turbulence power spectral contours obtained with this method are shown in Figure 
[TJ for the two diffusion tensors [eqs. (jlj) and ([5])] at two different times. Figure [2] shows 
the corresponding steady-state spectral contours (left panels) and power spectra at several 
propagation angles (right panels). 

As expected, the turbulence cascade in the parallel direction is initially suppressed. The 
Alfven wave period at the injection scale is given by Ta = 27v/k\\ t0 VA = 3.9 x lO 8 ^" 1 . It 
takes less than 10 Ta for the turbulence to reach the small scale boundary of the simulation 
domain. Steady-states are reached in about 10 Ta- This result is in conflict with simulations 
by Cho et al. (2002), where the turbulence power spectrum at a given k± is found cutting 
off exponentially with the increase of ku in the steady-state (Cho et al. 2002). It is not 
obvious to us what prevents the generation of waves in these simulations. Theoretically, 
if the turbulence cascades toward small scales following the critical balance proposed by 
Goldreich & Sridhar (1995), as confirmed by these simulations, it should be able to reverse 
cascade along lines (surfaces) of constant perpendicular wavenumbers through three wave- 
interactions. Our time-dependent numerical results in Figured] mimic such a reverse cascade 
process. Over a period of less than Ta (from the left panels to the right ones), high frequency 
nearly parallel propagating waves are generated by small scale fluctuations of the (quasi- 



7 When damping is included (see § 5 below), it dominates the cascade at high values of k (or small scales). 
We then choose an upper boundary for k large enough such that all injected wave energy is damped within 
the simulation domain. 
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Fig. 1. — Contours of the power spectrum W(k\\, k±, t) of the MHD range Alfven turbulence 
without damping at t = 3.46 x 10 9 Q- 1 (top left), 3.70 x 10 9 Q- 1 (top right), 1.16 x 10 9 Qp 1 
(bottom left), and 1.26 x 10 9 fi^ 1 (bottom right). The energy is injected uniformly at the 
surface of fe|| j0 = k±^ = 1-6 x 10~ 8 Q p /va with an injection rate F = 9.55 x 10~ w v A Q p . 
The logarithm to the base 10 of the contour levels in units of v A Q p 3 are indicated in the 
figures. The upper and lower panels show the results of the isotropic (-D||,|| = D± t ±) and 
anisotropic (-D||,|| ^ D±,±) diffusion models, respectively. The dashed lines indicate the 
critical balance relation discussed in the text. The power spect ra in the left panels are 



similar to those uncovered by MHD simulations (jCho et al. II2002I ). The reverse cascade in 



the perpendicular direction at high values of ku is evident in the right panels. 



2D) perpendicular component. This phenomenon is caused by the increase of the diffusion 
coefficient with the increase of the wavenumber and the suppression of the diffusion toward 
small parallel scales. At large scales, the wave energy can also diffuse slowly toward high 
values of ku. 



3.1. Results from Isotropic Diffusion Tensor 

Let us first consider the locally isotropic diffusion tensor (top panels). Note that the 
magnetic Mach number Ma = v(ko)/vA = 0.30, which means that our numerical solution is 
in the wave effect dominant domain. The steady-state power spectrum is nearly isotropic and 
is very similar to the IK spectrum between 10~ 6 Q p /va and 10~ 3 Q p /va, VV = Wo(k/k )~ 7 ^ 2+a , 



Fig. 2. — Left: Same as Figure [H but for the steady-state power spectral contours. Right: 
The scaled steady-state power spectra in several propagation directions indicated in the leg- 
ends. The Alfven Mach number Ma = 0.30 for the top panels and Ma = 0.23 for the bottom 
panels. The thick lines in the right panels are simulation results. The isotropic diffusion coef- 
ficient produces a nearly isotropic spectrum, which is very similar to the isotropic Kraichnan 
spectrum (thin horizontal lines), and the anisotropic diffusion coefficient produces spectra 
steeper than the Kolmogorov spectrum in directions other than the perpendicular direction. 
The thin straight lines (power-laws) are the approximate spectra based on equation (fT7|) . 
The spectral features at low values of k may depend on the geometry of the surface, where 
energy is injected into the system. 

where a(k) <C 1 is a slowly varying function of wavenumber k (a(k — > oo) — > 0). Substitution 
of this spectrum and Dy given in equation (jlj) into equation (jSJ) gives, the radial component 
of the slightly anisotropic energy flux density as 



where /i = cos 6 and the prime indicates derivative with respect to k. The total energy flux 
through spheres centred at the origin of the k space is then 



f{k) = (1/4tt)(7/2 -a- ka' \n(k/k ))W 2 k 5 [v k + v A fA 



-i 



8 



F 



(7/2 -a-ka' \n(k / k ))W 2 k 7 [ dfi[(v k + v A fi] 



-i 



Jo 



(7/2 -a- ka'\Ti(k/k ))W 2 k 7 vl 1 ln[l + v A /v k ], 



(9) 
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which should be independent of k in the steady-state case@ 

For a nearly isotropic solution, the constancy of F can be used to derive the weak 
dependence of a on fc: 

7/2 - a(k) - ka'\n(k/k ) = F v A /[W 2 k 7 ln(l + v A /v k )} (10) 

and the radial energy flux density becomes 

/= ^ , (11) 

4irk 2 (vk/vA + cosO) ln(l + v^/vy) ' 

which deviates from the k~ 2 dependence expected for exactly isotropic solutions for Ma < 1, 
especially near the perpendicular directions. So there is no exactly isotropic solution for this 
diffusion tensorj^ The steady-state solution we get here is just approx imately isotropic, and 



there is no feature associated with the critical balance proposed by iGoldreic h fc Sridhar 



(119951 ). It is similar to the steady-state fast mode solution given by IChandran I (120051 ) for 
coupled weak Alfven-fast mode turbulence. As expected, for v A , f — Fo/4iik 2 , the 

turbulence is isotropic. 

On the other hand, in the wave effect dominant regime if we assume an isotropic power- 
law spectrum with a = 0, the Alfven wave crossing process starts to affect the wave cascade 
at ku oc fc 3//4 , which is similar to the critical balance relation described by Goldreich & 
Shridhar (1995) and explains the fcn cx fc 3//4 scaling of the cascading wavevector front in the 
time-dependent power spectrum (Fig. [1]). This cascading wavevector front extends until it 



8 Note that the IK spectrum (VV cx k 1 1" 1 ) does approximately satisfy this constraint for va 2> v, and for 

v > va, 

F Q ~ (7/2 -a-ka' ln(fc//c ))W 3 / 2 fc 11 / 2 , 

one recovers the Kolmogorov spectrum. 

9 For locally isotropic diffusion tensors with their dependence on the wavenumber and the wave propagation 
direction sepa rable, one can al ways find an isotropic solution for the steady-state power spectrum in the 
incrtial range (jChandran 1 12005). In this case, the diffusion equation can be simplified as 



4nk 2 D(k, k, = F (k/k) , 



(12) 



where Fq gives the steady-state radial energy flux and is independent of k in the inertial range. Let FoQt/k) 
have the same propagation direction dependence as D, the corresponding formal solution is then given by 



W(k) = J dfcFo(k/fc)fc- 2 [47rL>(k/fc,/c,W)]- 1 



(13) 
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reaches the numerical boundary or the physical limit set by damping or non-MHD effects at 
large k. Because our eddy turnover timescale has different k and W dependence's than that 
given by Goldreich & Shridhar (1995), we obtain a different critical balance relation. For 



W oc k~ 11/3 with vl = Wk 3 ± and r NL = (vkk±) \ one can recover the Goldreich & Shridhar 
critical balance relation kn oc k 2 / 3 . 



The critical propagation direction, which separates the linear and nonlinear effect dom- 
inant regimes, satisfies cos6> c = v k /v A oc (k/ko)" 1 ^. Although the cascading energy flux 
is highly concentrated in the perpendicular directions with small values of ku, the diffusion 
process eventually carries wave energy to large values of k\\. The turbulence fills all possible 
wavevector space in the steady-state. The fraction of energy, which becomes linear Alfven 
waves, is given by the wave efficiency 

Vw= I fdtJ./[ fd f i = l-\n2/\n[l + v A /v k ], (14) 

J cos 8 C JO 

which approaches 1 as k — > oo. So the wave like component dominates at small scales as 
expected from the scalings of the cascade and wave crossing timescales with k mentioned 
above. 

When t>k > ^A, Vw < and the above formula for the wave efficiency is invalid. The 
turbulence then is dominated by nonlinear effects. One therefore can define a transition 
wave number k tT so that v(k tI ) = \W{k tr )kl v ] l l 2 = v A - Then 

MfcAt*) - (^) 2/3 ^ 11/3 (A:/A: tr )- 7 / 2+ ^< 1 (J^ * (k/k tt y 7/2+a , (15) 

where 

a(k) + ka'\n(k/k ) ~ 7/2 - (ll/3)(k tr /k) 2a / ln[l + v A /v k ] (16) 

and k tT — 3F /llv A . We note that a — > for k ^> k tr or v k <^ v A (IK spectrum), and a — > 
— 1/6 for k k tv or t\ ^> (Kolmogorov spectrum). The transition scale is proportional 
to the cube of the Alfven velocity and inversely proportional to the energy injection rate. It 
gives the coherent length of the magnetic field (Cho & Vishniac 2000) and may determine 
the maximum energy that charged particles can reach through resonant interactions with 
the turbulent magnetic field (Liu et al. 2008). The top right panel of Figure [2] compares the 
numerical results with the IK spectrum, W(k,F 0} v A ) ^ 0.22^ (11/3F ) 3 (k/k tr )- 7 ^ 2 (thin 
solid line). The deviations of the spectra in different wave propagation directions from the 
IK spectrum at k < 10~~ 6 Q p /v A shown in Figure [2] are likely caused by the cylinder geometry 
of the surface with fc|| = k± = k = 1.6 x 10~ 8 Q p /va, where turbulence energy is injected 
into the system. For the parameters of this numerical calculation, k tr = 2.60 x 10~ 10 Q p /va, 
and equation ffl~5l) with a = predicts an Ma = (fc tr //co) 1//4 — 0.36 at k = 1.6 x 10~ 8 f2 p /t>A, 
which is in agreement with the numerical results. 
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3.2. Results from Anisotropic Diffusion Tensor 



For the locally anisotropic diffusion tensor (jSJ) the diffusion coefficient in the perpendic- 
ular direction is higher than that of the isotropic diffusion model. Thus, for a given energy 
injection rate, in the steady-state case, the turbulence contains slightly less energy than that 
of the isotropic model. We find that for Ma = 0.23, the steady-state power spectral contours 
(Fig. [2j bottom panels) can be fitted with 



W(k) ~ 0.62 



3^0 
11 



2/3 



-11/3 



,tr 



+ 



k i 



,tr 



with a = 1, a± = 11/3, an = 4, ku tr = 0-5k± t t r , and k 



Hr 



2 \l/2 
tr) 



(17) 



3F /lk; A . The 



bottom right panel of Figure [2] compares this analytical fit (thin lines) to the numerical results 
(thick lines). The significant differences between the fit and numerical results at low values of 
k may also be related to the assumed cylindrical geometry of the injection surface. As evident 
for 9 ~ 7r/2 we have a spectrum very similar to the Kolmogorov spectrum, but away from 
the transverse direction the spectrum at high fc's becomes oc k~ A very quickly. In the limit of 
k > k tI , the angle averaged spectrum W(k) = fW(k)k 2 dtt oc k 1 ~ aa±+a±/a w = &r 7 / 4 . These 
results are similar to Chandran's steady-state solutions, but there are significant quantitative 
differences. The corresponding critical balance vxku = v^k leads to 



= 0.79kH 3 kl 11/e k 5 / 2 [l + 2 11 / 3 (k ]l /k ± )\k ± /k lltr ) 1 / 3 ]- 1 / 2 



0.79k%*k 



1/3,2/3 



This i s almost identical to the critical balance relation proposed by iGoldreich Sz Sridhar 
( 19951 ) 1^*1 The spectrum of the nonlinear effect dominant 2D component W^2d(^_l) = 27r/c^ J 



kj 2 , which is different from the result of IGoldreich fc Sridhar I (119951 ). As for the isotropic 
case, this is due to the difference in the scaling of the eddy turnover timescale with k and W. 



In our model v 2 oc Wk 3 , 



while v 2 oc ~Wk\\k 2 ± in Goldreich & Sridhar model. The fact that the 



spectrum of the 2D component is softer than the overall turbulence spectrum is also consis- 
tent with the assessment that the wave component dominates at small spatial scales. The 
energy flux of the 2D component F 2 D(k±) = 27rA:_|_ J v ^ k <vk D± t ±[dW(\si)/dk±]dk\\ oc kj 1 ^ 3 , 
which decreases with k±, and as the turbulence cascades towards small spatial scales, there is 
continuous energy conversion from the 2D perpendicular component to the wave component. 

On the other hand, based on the anisotropy of the steady-state spectrum, one may define 
a perpendicular component of the turbulence with k±/k± :tr > ku/ku tr , the corresponding 



spectrum W±(k±) = 2nk±_ J k 



W(k)d/c|| oc k ± . Due to the anisotropy of the dif- 



fusion tensor, this component develops quickly as can be seen from the time-dependent results 



10 



The critical balance only exists in the regime of <C va, which implies fcii <C k ~ fcj^. 
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(Fig. [T]). The perpendicular energy flux F±( k ±) = 27rk± J k k >k k D± : ±[dW(k) /dk±]dk\\ 

— 1/3 

oc k ± , which also carries most of the energy flux from large scales to small scales. The 
corresponding parallel component with k±/k± ttr < ku/k\\ t t r is dominated by linear waves and 
has a much steeper spectrum W^/cn) = f k±k t <k k±t 2nk 1 yV(k)dk± cx k~ 2 . The parallel 
energy flux F^(k\\) = J k±k ^ <k k±t 4:7rk±D\\ ! \\[dW(k)/dk\\]dk± oc fcjj" 1 . Since both the per- 
pendicular and parallel energy fluxes decrease with the increase of k, the total energy is not 
conserved, which reflects the inaccuracy of the analytical fit to the numerical results (Fig. [2]). 
The energy is better conserved with a slightly higher value of a±. Instead of pursuing even 
better analytical fits to the steady-state spectrum, which most likely will be complicated, we 
focus on numerical results in what follows. 

The most notable difference between these results and recent theoretical and numerical 
studies of Aflven turbulence is the presence of a prominent parallel wave component in 
our model. Such a component is found to be absent or dynamical unimportant in many 
studies (Galtier 2006; Clio et al. 2002; Goldreich & Sridhar 1995). On the other hand, 
these parallel wave modes may play a dominant role in the selective acceleration of ions from 
the background plasma (Mason et al. 2002; Liu et al. 2004, 2006) and scattering of high 
energy particles. They have been the subject of extensive studies (Yan & Lazarian 2002; 
Farmer & Goldreich 2004; Lazarian & Beresnyak 2006). Since the diffusion process tends to 
fill the wave vector space with energy, the presence of a strong parallel wave component is 
not a surprise in our model. They are also not subject to the turbulence damping studied 
by Farmer & Goldreich (2004) and Lazarian & Beresnyak (2006). This damping is caused 
by the interactions of the waves with the quasi-2D strong turbulence component in the 
perpendicular direction, which makes the waves cascade toward small perpendicular scales. 
Our parallel wave component is produced by the cascade (or reverse cascade) itself. So the 
cascade process generates these waves instead of damping them. With the strong turbulence 
theory proposed by Goldreich & Sridhar (1995), high values of k\\ can be reached through 
the critically balanced cascade k\\ oc k 2 / 3 . The turbulence power spectrum is given by 
W oc g{k\\/h±). If the turbulence can cascade along lines of constant values of k± as 
required by the three-wave resonance conditions, one would expect that at a given value of 
fey, the spectrum is flat with the increase of k± until the critical balance condition is satisfied. 
Therefore W oc kj 10 ^ 3 oc k^ 5 . This result is similar to the result of our anisotropic cascade 
model. However, they are quite different from the power spectrum derived from MHD 
simulations by Cho et al. (2002), where it is found that the spectra cut off exponentially 
toward small parallel scales. It remains to be seen what prevents the presence of the parallel 
wave component in these simulations. 
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4. Turbulence Cascade beyond the MHD Regime 

So far we have assumed the simple MHD dispersion relation uj = kvA cos 9, which is valid 
on large scales. For a more exact treatment of both cascade and damping we need to go be- 
yond the MHD regime and use a more complete description of the dispersion relation, which 
is well known. Linear waves in collisionless magnetized plasmas have been studied extensively 
with the linear Vlasov theory (Andre 1985). Here, for the purpose of comparing with the ob- 
served turbulence in the solar wind, we focus on the dominant Alfven-cyclotron fluctuations. 
Observations suggest that the Alfven turbul ence can cascade to spatial scales below the MHD 



regime and the kinetic effects are important (IDenskat et al. Ill983l ; iLeamon et al. 1119981 11999 



Bale et al. 1 120051 ; lOsman fc Horbury 1 120071 ) . For fully ionized plasmas with the solar abun- 
dance, the dispersion relation deviates from the simple MHD relation significantly near the 
4 He nucleus (or a-particle) cyclotron frequency Q a = 0.5Q P = 0.076(27r)(i?/10nT)Hz (the 
Alfven-cyclotron fluctuations exist only below this frequency). 

Our diffusion tensor is constructed from the wave dispersion relation w r (k), where the 
subscript "r" refers to the real part of the frequency. (As discussed in §2.2 and below, the 
imaginary part determines the damping rate.). For an accurate determination of u r we use 
the WHAMP code to obtain the "hot-plasma" Alfven wave dispersion surface (Ronnmark 
1982). The thin lines in Figure [3] show this relation for different angles of propagation under 
typical solar wind conditions (temperature T ~ 4 x 10 4 K). These match very closely the 
"cold plasma" (T = 0) dispersion relation, shown by the thick lines (see Appendix A), 
up to k ~ O.Q^I p /va with less than 30% deviation at highest frequencies!"] At a higher 
(lower) values of temperaure the deviation will start at a lower (higher) value of k. In 
what follows we use the cold plasma dispersion relation without loss of much accuracy but 
with a great gain in computational simplification because these dispersion relations can be 
described analytically. This greatly simplifies the construction of the diffusion tensor at all 
the wavevector grid points. 

For Alfven waves, the turbulence cascades preferentially in the perpendicular direction. 
This is still true when we consider the kinetic effects on the plasma wave dispersion relations. 
However, as mentioned above in the MHD regime the Alfven waves have a well defined 
wave crossing time. The kinetic effects introduce ambiguities on the definition of rw- One 
therefore can identify r^ 1 with either ov(k) or v gr -k. Because the Alfven wave packet becomes 



n Note that the hot plasma dispersion curves terminate at some finite values of k. This is because currently 
there are no simple schemes, which can lead to reliable dispersion relation (and damping rate) in the high 
frequency range, where damping dominates the wave effect. At such high frequencies thermal damping 
dominates and, as we will see below, the turbulence spectrum cuts off sharply. 



-22 - 




Fig. 3. — Dispersion relations (the real part of the wave frequency) of Alfven-cyclotron 
fluctuations in several wave propagation directions. The thick lines give the cold plasma 
relations and the thin lines are obtained numerically with the WHAMP code. We assume 
Q p = 0.096(27r) Hz (B = 6.3nT), the electron and ions temperatures k B T e = k^T p = k B T a = 
3.35 eV, and densities n e = 12 cm -3 (/3 P = 0.41), n p = 10 cm -3 , n a — 1 cm -3 , which are 
typical for the solar wind (Leamon et al. 1998, 1999). Note that the Lauren series are used 
in the code to approximate the plasma dispersion function Z, that is valid when the damping 
rate, i.e., the imaginary part of the wave frequency Ui = T, is much less than the real wave 
frequency u r . Thus, the solution at high frequencies may not be valid. The missing segments 
of the dispersion relation in the figure are related to this caveat. 



stationary near Q a (where k — > oo), one may favor the latter choice that implies negligible 
wave propagation effects. It is not obvious which one of these two choices is physically 
more reasonable. We first try both forms still without the damping (or leakage) term. The 
corresponding steady-state spectra are shown in Figure HI The spectra are nearly identical 
at large scales (the spectra for r^ 1 = uj are shifted upward by a factor of 5 for the illustrative 
purpose) and, as expected, both contain a spectral break at kuv^ ~ where the wave 
dispersion surface starts to deviate from the MHD relation. These breaks, though appealing, 
only have a spectral inde x change of less t han one. They may explain a few observations 
of solar wind turbulence (IBale et al. II2005I ). However, neither of the se breaks can account 



for the more commonly observed steep spectra at high frequencies (IDenskat et al. 1 11983 



Leamon et al. 1 1 19981 ; ICranmer fc van Ballegooijen 1 120031 ; lOsman fc Horbury 



20071) • In the 
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kinetic effect dominant regime, the wave propagation effects vanish for r w x 
spectra approach, once again, the Kolmogorov spectrum at high values of k. 



v gr -k, the 



We will show below (§|5j) these two cases cannot be distinguished with most observations 
of the solar wind turbulence since the strong thermal damping cuts off the spectra at the 
wavevectors before the spectral deviations due to the kinetic effects shown above set in. 
[For other conditions, e.g., with lower temperatures, the break due t o kine tic effects may 
appear before the damping cuts off the spectra, see e.g., iBale et al. I (120051 ).] The kinetic 
effects, however, play crucial roles in the thermal damping processes. In what follows, we 
leave the subtle differences between these two possible choices of the wave coupling times to 
further study, and adopt the first choice, i.e., Tw = obtained from the exact cold plasma 
dispersion relation for both the isotropic and anisotropic diffusion tensors [eqs. PJ and (j5j)]. 
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Fig. 4. — The steady-state wave power spectra in several wave propagation directions. The 
energy is injected at k\\ = k ±0 = 2.01 x 10- 4 a p /v A with W(k ) = 1.44 x l0 9 v A tt p 3 . The 
corresponding Mach number M A = 0.18. The upper and lower lines are for r-^ 1 = to r (k) and 
k • v gr , respectively. To avoid confusion the former has been shifted upward by a factor of 5. 
The spectral breaks at k\\V\ — 0.5Q P in both cases are due to the kinetic effects on the wave 
dispersion relation. 
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5. Turbulence Spectra in the Dissipation Range 



It is well-known that the d amping of waves in a m agnetized plasma by thermal back- 
ground particles is not isotropic (I Gary Sz Borovsky l l2004h . In combination with the anisotropy 



of the turbulence cascade, we would expect strong anisotropics in the turbulence spectrum 
in the dissipation range. These anisotropics have significant effects on the energizing of 
background particles by turbulent plasma waves and have broad implications on the energy 
dissipation in collision-less astrophysical plasmas. Although these issues have been well- 



Cranmer & van Ballegooijen 


2003; 


Howes et al. 


2007. 


2008a) 



Compared with previous studies, our model has less assumptions with most of the related 
physical processes treated more self-consistently. We present the corresponding results in 
this section. 



5.1. Start of the Dissipation 

As mentioned above the damping rate can be calculated from the WHAMP code for 
a thermal plasma with T(k) = Wj(k), the imaginary part of the wave frequency. To our 
knowledge there are no equations or programs for damping by an a rbitra ry nonthermal 



particle distributions except the recent result from iPetrosian et al. I (120061 ) for the TTD 
damping of fast mode waves in low beta plasmas. As shown in this paper the damping rate 
increases with increasing k more rapidly than the cascade rate so that it becomes important 
beyond some critical wavevector k c obtained from T(k c ) = r c ^(k c ). And one expects sharp 
cutoff of the turbulence spectrum for k > |k c | at a given wave propagation direction, which 
corresponds to the dissipation range. In what follows we first follow this procedure assuming 
a thermal plasma to estimate the dissipation range turbulence spectrum and then show exact 
numerical results from the solution of the wave equation including both the diffusion and 
damping terms (with the leakage term still ignored). 

For the locally isotropic diffusion model, the cascade rate is uniquely defined as t~J. = 
t nl/(^ + t nl/ t w)- The left panel of Figure \5\ shows this cascade rate (from the above 
numerical results with damping excluded) and the thermal damping rate (for the plasma 
in Figure [3]) in several wave propagation directions. We see that with the increase of k 
the damping rate increases much faster than the cascade rate. The cutoff wave number k c 
is given by the intersection of these two rates. The solid line in the middle panel shows 
the dependence of k c on 6. The dotted and dashed lines give contours of constant cascade 
and damping rates, respectively. These anisotropic damping and cascade make k c change 
by more than one order of magnitude with the change of the wave propagation direction. 
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This is distinct from the ID turbulence model, where the dissipation range covers a very 
narrow spatial scales with the Reynolds number on the order of 1. The right panel shows 
the dependence of k c on (3 P and the turbulence intensity. Higher values of (3 P lead to stronger 
damping and lower values of k c . With the increase of the turbulence intensity, the dissipation 
range shifts to smaller spatial scales. For the model with the anisotropic diffusion tensor, the 
cascade rate in the parallel direction is different from that in the perpendicular direction. 
The model expectations are also more involved. However, numerical evaluations of both 
models are straightforward. 




Fig. 5. — Left: The wave cascade rates obtained without damping (thick) and damping 
rates thin) as functions of k for three wave propagation directions. The properties of the 
plasma are the same as in Figure [3j The energy is injected at k\\o = k± = 2.01 x 10~ 4 fip/tA 
with F = 1.89 x 10~ 3 t> A fi p , which implies W(k ) = 6.62 x 10 10 v A ^V 3 - The corresponding 
Ma = 1.23. Middle: The dotted and dashed lines show a sample contour of the above cascade 
and damping rates in the wavevector space (k, 6), respectively, for r _1 = 10 _4 f2 p . The critical 
wavenumber k c defined as the location where these two rates are equal is indicated by the 
solid line. At a given 9, u>i increases monotonically with k. The turbulence spectrum should 
cut off at k c , which varies by more than one order of magnitude with the change of the wave 
propagation direction. The nearly parallel propagating waves are damped through cyclotron 
resonances at uo ~ Q a - Obliquely propagating waves are subject to the TTD and LD. The 
nonlinear effects dominate in the nearly perpendicular directions. Right: The dependence 
of k c on j3 p and Ma- Solid, dashed, and dotted lines show the cases with (Ma = 1.23, 
/3 P = 0.41), (M A = 0.77, (3 P = 0.41), and (M A = 0.77, f3 p = 1.0), respectively. With the 
increase of the turbulence intensity, the cascade rate increases leading to higher values of k c . 
The damping rate increases with (5 P so that a higher j3 p leads to lower values of k c . 
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5.2. More Exact Spectra 

We now present some numerical results on the spectrum and anisotropy of the turbulence 
in the dissipation range where we include the thermal damping term using the WHAMP code 
as described in § 12.21 The left and middle panels of Figure [6] show the steady-state turbulence 
power spectrum for the isotropic diffusion model. For the plasma parameters chosen here 
the break due to damping comes before the spectral break associated with the kinetic effects 
shown in the previous section. We there fore do not expect sp ectral features caused by the 



kinetic effects on the turbulence cascade (IStawicki et al. 1120011 ). The spectra cut off sharply 



when the damping dominates, which is consistent with the results of the previous ID models 
(Miller et al. 1995, 1996; Li et al. 2001; Stawicki et al. 2001). This is mainly due to the 
non-linear nature of the diffusion equation. Once the damping rate is high enough to make 
the turbulence spectrum deviates significantly from the inertial range spectrum, the cascade, 
whose rate decreases with the decrease of the local turbulence power in the wavevector space, 
is suppressed, driving the cascade and damping processes out of balance quickly with the 
increase of k and resulting in a sharp spectral cutoff. Therefore, it is very difficult to produce 
a broken-power law spectrum or even a gradual cutoff through the thermal damping with 
the ID model except in some unusual scenarios, where the damping processes themselves 
are al so nonlinear and the da mping rate scales the same way with k and W as the cascade 



rate (IMarkovskii et al. 1120061 ). We also don't see the gradual spectral cutoff due to th e 



gradual onset of damping with the increase of k as suggested by iGary fc Borovsky I (12004 ). 



Such more gradual c utoffs may appear with different properties of the background plasma 



flHowes et al. Il2008ah 



As expected, the actual location of the cutoff wavevector depends on the wave propaga- 
tion angle. As a consequence of weaker damping and stronger cascade in the perpendicular 
directions, the perpendicular wave spectra cut off at scales (wavenumbers) about one order of 
magnitude smaller (larger) than that of the parallel waves. Using the equality of the cascade 
and damping rates as the criterion for spectral break one would expect strongest damping for 
waves with 9 ~ 10°; the minimum of k c vs. 9 plot (in middle panel of Fig. 5). However, the 
valley here is smoothed out due to the nature of the diffusion processes. Waves with higher 
and lower values of 9 can cascade to smaller scales and induce waves with 9 ~ 10° through 
the diffusion processes. The isotropy of the diffusion tensor makes the spectrum of waves 
subjected to heavy damping cut off at a value higher than k c . Note also that the cascade 
rate in Figure \5\ is evaluated without the damping. The actual cutoff wavenumber should 
be lower than k c in general. The right panel shows the spectral contours for the anisotropic 
diffusion tensor. The parallel component is damped even more quickly because the turbu- 
lence spectrum is dominated by the perpendicular component whose damping sets in within 
a relatively narrow wavenumber range. The turbulence is therefore strongly anisotropic in 
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Fig. 6. — Left: Same as Figured with the thermal damping included. The turbulence injec- 
tion and properties of the background plasma are the same as that in the left panel of Figure 
|5j As expected, the power spectrum cuts off sharply at high values of k, and the location 
of the cutoff wavenumber k c is very sensitive to the wave propagation direction. Middle: 
The corresponding contours of log 10 (W / /c 11//3 ) (solid lines, starting at -2.5 and decreasing in 
steps of 0.5) of the Left panel in k — 9 plane. The spectrum breaks near —2.5. The dashed 
line indicates k c obtained from the prescription given in Figure [51 The diffusion process 
smoothes out the valley near 9 = 10° in k c (9). Right: Same as the Middle panel but for the 
anisotropic diffusion tensor. The turbulence is dominated by the perpendicular component, 
which is damped over a narrow wavenumber range. 



the dissipation range. 



5.3. Total Spectra 

The angle averaged turbulence spectra are shown in the left (as a function of the 
wavenumber) and middle (as a function of the frequency) panels of Figure [7J The spec- 
tra of the isotropic diffusion model mimic broken-power laws with softer dissipation range 
spectra. The dissipation range extends nearly one order of magnitude in both wavenumber 
and frequency, which may explain the observed broken power law spectrum of the solar 
wind Alfven turbulence. However, due to the low value of fl a , the break frequency is about 
one order of magnitude lower than the observed value. Although the anisotropic diffusion 
model has a broader dissipation range as shown in the right panel of Figure El due to the 
dominance of the perpendicular component the angle averaged spectrum can be fitted with a 
single power law, which cuts off when the damping of the perpendicular component becomes 
significant. The spectrum of the anisotropic diffusion model therefore appears to cut off 
much sharply with a narrow dissipation range, which is similar to that of the ID diffusion 
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Fig. 7. — Left Angle averaged turbulence power spectra as functions of k for the spectra 
in Figure [6j The isotropic diffusion tensor (solid) gives a gradual cutoff with the dissipation 
range covering about one order of magnitude in the wavenumber as expected from the Middle 
panel of Figure O The overall spectrum mimics a broken power law with a spectral break 
near k ~ 0.5fi p /t>A (for the assumed plasma parameters). The anisotropic tensor (dashed) 
gives a much sharper cutoff at k ~ l/p p m agreement with the results of Howes et al. (2007), 
where p p is the gyroradius of the thermal protons. Note that v\ ~ t> t h. p for the assumed 
plasma parameters. Middle: Same as the Left but plotted as a function of the frequency. 
Both spectra cut off at Q a though the cutoff of the anisotropic model is sharper. The 
isotropic model has a broader dissipation range with the overall spectrum mimics a broken 
power-law with a sharp high frequency cutoff. Right Same as the Middle panel but with the 
Doppler shift due to the relative motion between the solar wind and the spacecraft included. 
The angle between the magnetic field and solar wind velocity 6bv = 38° and the solar wind 
speed Vsw — 517km/s. As in the Left panel, the spectrum of the isotropic diffusion tensor 
can be fitted with a broken power law (dotted) with the break frequency v^f = 0.2Hz and 
the lower and higher frequency spectral indexes 71 = —1.7 and 72 = —3.1, respectively. 



models. The spectrum of the isotropic diffusion model is also slightly harder and higher than 
that of the anisotropic model in the "inertial" range, which is in agreement with our study 
of the Alfven turbulence in § [31 

In summary, the above results demonstrate the possibility of a rich variety of spectra for 
turbulence determined primarily by the plasma parameters; the density n , temperature T, 
magnetic field B, source size L, and the inje cted energy flux Fq and wavenumber ko ~ L^ 1 
(ISmith et al. 1 120061 ; iMarkovskii et al. 1 120081 ) . These parameters can be derived from MHD 
simulations, which, in combination with the present diffusion model, can be used to study the 
free energy dissipation processes in magnetized plasmas quantitatively. The uncertain part 
of this procedure is the choice of the diffusion tensor. As discussed above we have constructed 
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tensors that produce many features of turbulence spectra obtained with simulations based 
on more basic physics. Such simulations are too cumbersome or limited to be useful for 
comparison with observations of astrophysical sources. But they are important in guiding us 
to derive the appropriate form of the diffusion tensor that lies at the core of our approach. For 
a given diffusion tensor, the above discussion demonstrates that one can readily (that means 
no need of supercomputing power) produce spectra for turbulence, for varied background 
plasma conditions, which can then be used to determine the heating of the background 
plasma and acceleration of particles and the resultant radiations. The validity of the chosen 
diffusion tensor form can be esta blished by comparing the model predictions with existing 
observations or some simulations (IHowes et al. 1120081 ; iGary. Saito. &: Li II2008I ). In the next 
section we compare the predictions of our models with the observed spectra of magnetic 
fluctuations in the solar wind. 



6. Application to the Solar Wind Turbulence 

To apply the model to observations of the solar wind turbulence, one must take into 
account the Doppler shift due to the relative motion of the spacecraft with respect to the 
solar wind (Leamon et al. 1999). The observed frequency spectrum P(v) is related to the 
turbulence spectrum as: 

P{u) = j W(k)5 ji- [k • V sw + u(k)] - v\ dk (19) 

where v is the spacecraft-frame frequency and S is the Dirac-5 function. The right panel of 
Figure [7J shows the P(y) for the angle between the magnetic field and solar wind velocity 
6bv = 38° and the solar wind speed Vsw — 517 km/s. The rest of the model parameters 
are the same as the middle panel, which would be the observed spectrum for Vsw — 
km/s. [The spectrum with respect to the wavenumber k in the left panel can be obtained 
by replacing the argument of the 5 function with (k — k'), where k! is the variable of the 
integration.] 

All these spectra can be fitted with broken power-laws with steeper high end declines. 
To understand this result, we may make the approximation: W(k) oc k~ n ^ 3 Q[k c (9) — k\. 
Then the spectrum on the left panel of Figure [7] is given by P(k) oc k~ 5 ^ 3 cos9 c (k), where 
9 c (k) is the cutoff propagation angle in the dissipation range and is given by the inverse of 
the function k c (9) shown in Figure which increases monotonically with 9, except for a 
small range near 9 = 0. Clearly, P(k) in the dissipation range is determined by k c (9). Since 
k c (9) covers more than one order of magnitude, the same is the dissipation range of P(k). 
The spectrum in the middle panel (Vsw = 0) is obtained by changing the variable from k to 
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frequency according to the dispersion relation o;(k). In the MHD regime u;(k) = VAku, then 
we have P{u) oc z/~ 5//3 cos 8//3 9 c (k). Beyond the MHD regime u;(k) flattens and approaches 
Q a as k — > oo. As a result -P(z^) cuts off at Q a /2ir. For the spectrum in the right panel we 
must carry out the integration in equation f[T9~j) . For solar wind conditions in general, and 
for the observations we compare with below in particular, the Alfven velocity (va = 36.8 
km/s) is more than one order of magnitude lower than solar wind velocity (Vsw — 517 
km/s). One may then ignore the w(k) term in the argument of the ^-function in the above 
integration. Then for an isotropic power-law turbulence spectrum, the shape of P(v) will 
be identical to P(k) in the inertial range. In general, P(v) oc z/~ 5 / 3 J CO 
9)dcos9, where k = 2irv/Vsw cos(9bv — 9 C ). The dissipation range of P{y) starts at v d ~ 
Vswk c (0°) cos6 , sy/27r ~ 0.1 Hz for the parameters in Figured where we have assumed that 
k c (9) increases monotonically with sin#. 

The break frequency of the broken power law fit will be at a higher frequency; v b f ~ 
Vswk c (9 W) /27T '. By fitting P(v) with a broken power law model down to the power spectral 
level of 2 x 10 -4 nT 2 /Hz, we obtain the dotted line in the right panel of Figure [7] with 
Ubf = 0.2 Hz, and the lower and higher frequency spectral indexes 71 = —1.7 and 72 = —3.1, 
respectively. (For clarity the dotted line is shifted upward by a factor 2.) This broken power 
law spectrum is very similar to those observed in the solar wind (Leamon et al. 1998, 1999). 

6.1. Fit to Observed Spectra 

As an example, we fit the observed solar wind turbulence spectrum at 2200UT, January 
11, 1997 at 1 AU from the Sun (Leamon et al. 1999), as shown by the solid line in Figure 
El The observation shows that 9 BV = 38°, V sw = 517 km/s, B = 63^G [Q p = 0.096(2vr)Hz], 
and /3 P = 0.48. The Alfven velocity (or the plasma density or temperature) is not given 
by the authors. By adjusting v\ and the energy flux Fq carried by the turbulence, we 
obtain the best fit to the observed spectrum for the locally isotropic diffusion tensor model 
as shown by the dashed line in the Figure El Leamon et. al. (1999) fit the observation with 
a broken power law model with indexes 7; = —1.67, 7^ = —2.91, and a break frequency 
Vb = 0.235Hz. Our simulated spectrum, when fitted by a broken power law model down to 
2 x 10- 4 nT 2 /Hz, gives 71 = -1.67, 72 = -2.97, and u bf = 0.200Hz. F is determined by 
the turbulence spectrum in the inertial range. With va as the only free model parameter, 
the model reproduces not only the observed break frequency v b f, but also the power- law 
index in the dissipation range. This is the most appealing success of the model. Our model 
fit gives v a = 34 km/s, n e = 14 cm -3 for k^T e = 3.35 eV, which can be compared with 
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Fig. 8. — Observed interplanetary turbulence power spectrum (solid line; extracted from 
Figure 1 of Leamon et. al. 1999) compared with our model spectrum (dashed line; scaled 
downward by a factor of 5 for clarity) based on the locally isotropic diffusion tensor and 
the model parameters B = 6.3nT, (3 V = 0.48, 9bv = 38°, Vsw — 517km/s given by the 
observations. The Alfven velocity is not well determined. Our best fit model has an Alfven 
velocity of 34 km/s. The observed spectral flattening at high frequencies is related to a 
spacecraft noise. As evident this model with essentially only one free parameter provides 
an excellent fit to the observed spectrum. We note that this model is incompatible with 



the kn own properties of Alfven turbulence, such as those studied by iGoldreich &: Sridhar 



(Il995[ ). 



observations^ 

We also find that the model based on the anisotropic diffusion tensor cannot fit the 
observed spectrum very well. This is primarily due to the dominance of the perpendicular 
component and the much sharper cutoff in the dissipation range as shown by the dashed lines 
in Figured for the angle averaged spectra. We also consider the Kolmogorov phenomenology 



12 The cascading constant C in equation (j4|) should be considered as another parameter. There is also an 
uncertainty in the definition of the eddy speed. For this simulation, we set C — 1. A complete measurement 
of the propertie s of solar wind plasma s and the turbulence carried by them can be used to determine C and 
test the model |Yeung fc Zhou Ill997h . 
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with the isotropic cascade rate r^J. oc fc 2//3 in the inertial range. The damping rate is very 
sensitive to k but has relatively weak dependence on 9 (Fig. [5]). So k c , which is identical to 
the contour of the damping rate, has a weaker dependence on 9 than the isotropic diffusion 
model studied above. As a result, the angle averaged spectrum has a sharper cutoff at high 
frequencies similar to the anisotropic diffusion model. 

Considering the kinetic effects on the wave cascade and dampin g in the perpendicular 



direct io n along the li n e defined by the critical balance proposed by iGoldreich fc Sridhar 



(1995 


), 


Howes et al. 


f 


2007, 


2008a) 



broken-power like spectrum with a cutoff at small spatial scales. In this model, the kinetic 
Alfvev wave (KAW) has a steeper spectrum than the Alfven waves. However, to have 
significant change in the spectral shape between the Alfven and KAW waves, the electron 
temperature needs to be much higher than the proton temperature to make the kinetic effects 
significant. This is quite different from our model, where a broken-power like spectrum is 
produced by the anisotropic damping and the inhomogeneous cascade in the wavevector 
space. The broken power-law spectrum they produced is with respect to the perpendicular 
component of the wavevector. To compare with solar wind observations, one needs to take 
into account the Doppler effects. The KAW starts at kj_p p rs 1. The corresponding Doppler- 
shifted frequency 

v*k.V sw /2n>k ±Pp ^- 2 ^. (20) 

For typical solar wind conditions with /3 P < 1 and V$w ~ 10fA, this frequency is more than 
10 times higher than the proton gyro-frequency, which disagrees with the observed fact that 
the break frequency is usually slightly higher than VL p /2k (Leamon et al. 1998, 1999; Bale 
et al. 2005). With the ID model, Stawicki et al. (2001) argues that the broken power-law 
spectrum of the solar wind turbulence may be attributed to the dispersive effects alone as 
the turbulence cascades from the Aflven wave domain to the whistler wave regime. However, 
the diffusion coefficient must be prescribed properly to reproduce the observed spectrum. 

From the above discussion we conclude that the model based on the locally isotropic 
diffusion tensor provides the most natural explanation of the observed solar wind turbulence 
spectrum. All other models need fine-tuning of the model parameters or relevant physical 
processes to give acceptable fits to the observed turbulence spectrum. 



6.2. Spectral Features and Model Parameters 



The previous section shows the success of the model in accounting for a specific observa- 
tion interval. To understand the nature of the observed broken power spectrum, statistical 
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studies of a large samp l e of events are necessary. Some of these stud ies ( ILeamon et al 



1998; 


Smith et al. 


2006; 


Markovskii et al. 


2008; 



Hamilton et al. II2OO8I ) have uncovered the 



nonlinear nature of the dominant dynamics and have shown that several factors may affect 
the spectral shape in the dissipation range. For a qualitative comparison with these ob- 
servations, in this section we explore the dependence of the spectral features on the model 
parameters. 
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Fig. 9. — Dependence of the angle averaged spectrum on 9 the angle between the magnetic 
field and solar wind velocity (the same as the 9bv m the text). Left: Two sample spectra 
with Ma = 1.23, (3 P = 0.41, and two values of 9bv as indicated in the legend. The other 
parameters are the same as in the right panel of Figure [7J Broken power law fits to these 
spectra down to 2 x 10~ 4 nT 2 Hz _1 give a break frequency Vbf = 0.17 Hz for 9bv — 15° and 
= 0.24Hz for 9bv — 60°. Right Dependence of and 72 on 9bv- Two lower limits of 
the power spectrum are considered while fitting the spectrum with broken power laws. The 
cross signs are obtained by fitting the model spectrum down to a turbulence power level of 
2 x 10 _4 nT 2 Hz _1 , and the circle signs correspond to a low limit of 2 x 10~ 7 P(z/ = 10 _3 Hz). 

As indicated above, most observed spectra can be fitted with three parameters; 71 
that is always nearly equal to -5/3, 72, and the break frequency v b f. We have explored 
the dependence of the latter two on the observable parameters using the locally isotropic 
diffusion tensor scheme. Figures [9j [101 and [TT1 show the dependence of the spectrum on ^bv, 
plasma beta /3 P , and Mach number Ma, respectively. The other parameters of the plasma 
and the injected turbulence are the same as that in the right panel of Figure [7J unless specified 
otherwise. To describe the spectrum quantitatively, we fit the numerically calculated spectra 
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to a broken power-law model down to a turbulence power level of 2 x 10~ 4 nT 2 Hz _1 and 
2 x 10 _7 P(z/ = 10~ 3 Hz). The former is appropriate to compare with observations made by 
instruments with a given sensitivity to the magnetic field turbulence power. The latter is 
appropriate for observations with a given dynamical range in the turbulence power spectrum. 
In general, and 72 are very similar for these two fits, and as stated above, 71 is always very 
close to —5/3 for the parameter space explored. The most distinct feature of these results is 
the correlation between the two observables u b j and 72 (Fig. [12]). A higher always comes 
along a softer dissipation range spectrum. The correlation caused by variations in /3 P and Ma 
is related to the kinetic effects of the waves beyond the MHD regime. The damping rate is 
very sensitive to (3 P , as can be seen from the right panel of Figured The change in Ma leads 
to different dissipation scales, where the kinetic effects are different. The resultant correlation 
between v h $ and 72 is mostly caused by the fact that the cascade rate becomes more isotropic 
near fl a (the left panel of Fig. [5]). With the decrease of the dissipation scales caused by 
either higher values of Ma or lower values of /3 P , the dissipation range covers a narrower 
spatial extension, giving rise to a softer dissipation range spectrum. The dependence on 
Obvi 011 the other hand, is mostly due to the Doppler shift. Since ~ Vsw^c^bv)/^ 71 , v bf 
increases with sin#BV> and the spectral range beyond v^j becomes narrower leading to lower 
values of 72. These explain the nearly identical correlation between 72 and Ubf caused by 
variations in [3 P and Ma while a quite distinct correlation due to changes in 9bv, as shown 
in Figure [121 
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Fig. 10. — Same as Figure[9]but for the dependence on /3 p . Here Ma = 1.23 and Obv = 38°. 
However, the dependence of the spectrum on 9bv in our model is much weaker than 
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any quasi-2D turbulence models with the energy cascading only in the perpendicular direc- 
tion (IHowes et al. 1 120071 . l2008at ). If most of the wave energy resides in the perpendicular 



direction, as suggested by the critically balanced cascade proposed by Goldreich & Shrid- 
har (1995), the k c contour will be highly prolonged in the perpendicular direction and 9bv 
will str ongly affect the break freq uency. The observations by Leamon et al (1998, Fig. 6) 
and by lOsman fc Horburv I (120071 ) appear to favor the isotropic diffusion model used here 
( iHamilton et al. II2008I ). The dependence of the spectrum on the solar wind speed is rela- 
tively simple. Obviously, a faster solar wind gives a higher break frequency. Clearly, more 
observations and/or systematic analyses of existing data are required to test the validity of 
these models. 
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Fig. 11. — Same as Figure[9]but for the dependence on Ma. Here /3 P = 0.41 and 6bv = 38°. 

In fact there are already some preliminary results showing some weak correlations (with 
lots of scat ters) between the sp ectral features and other parameters. For example, the 



Figure 3 of iLeamon et al. I (119981 ) shows a correlation between the break frequency and the 
proton cyclotron frequency. For a given thermal pressure, higher values of proton cyclotron 
frequency imply lower values of /3 P , and as shown in Figure [TD1 we expect higher break 
frequencies as observed. Their Figure 6 shows two events with distinct values of 9bv but 
similar values for the rest of the parameters. This observation suggests the dissipation range 
spectrum is harder for higher values of 9bv m contradiction with our Figure [9j However, 
their spectrum with a higher value of 9bv shows prominent noise at high frequencies so that 
the spectral break is obscured. Their Figure 2 suggests that the dissipation range spectrum 
is softer for higher values of proton temperature T p , which implies higher values of (3 P for a 
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given Alfven speed. While we expect the opposite, we notice that the higher values of f} p 
likely result from higher values of Ma, as indicated by the convergence of their inertial range 
spectral index with the increase of T p and simultaneous increas es of Ma and (3 P can make 
the dissipation range spectrum change in either direction. Indeed, ISmith et al. I (120061 ) found 
that the dissipation range spectrum becomes softer with the increase of the energy cascade 
rate in the inertial range. According to our model, the energy cascade rate is determined by 
Ma and k . Thus this observation has the same trend as the dependence of 72 on Ma shown 
in Figure [TTJ 
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Fig. 12. — The correlation between v^f and 72 caused by variations in 9 (crosses), f3 v (circles), 
and Ma (pluses). The model parameters are the same as those in the right panel of Figure 

m 

Leamon et al. (1999, Fig. 2) showed a weak correlation between fc^iss (the same as 
our k c ) and the inverse of the proton gyro-radius Rl (the same as our p p ). For damping 
dominated by cyclotron resonances, fcdissCs ~ Q p , where Cs is the sound speed. Therefore 
Ubf oc fcdiss oc fl p / Va^/Wv-i w hich is similar to the trend we find in Figure [TOl In general, higher 



values of Rl imply higher values of /? p and longer dissipation scales. ISmith et al. I (120011 ) 



studied a solar wind interval with very high Aflven speeds and found that the cutoff frequency 
decreases with the increase of va- They concluded from this that the dissipation is related 



13 The inertial range spectrum can be fully developed in shorter periods for turbulence with higher values 
of M A . 
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to the ion inertial length v\/fl p , instead of Rl- We note that the turbulence intensity was 
very low during this interva l , whic h implies low values of u^f according to our Figure [TT1 More 
recently. iMarkovskii et al. I (120081 ) correlated the break frequency u^f with several parameters 
and found some weak correlations. For example, their Figures 4 and 7 also show a decline 
of break frequency with the plasma beta; Vbf oc (3~ x with x ~ 0.24, in agreement with our 
Figure [TOl On the other hand, their Figure 6 can be interpreted as V\,f oc 1/Ma, which has 
the opposite trend to what we find in Figure [TTJ 

We only have qualitative comparisons of these observations with our model here since 
the turbulence spectrum depends on several parameters, Ma, 9bv, Vsw, Pp, etc. Simple 
correlations between 72, v h $ and any one of these quantities are not expected because varia- 
tions caused by the other parameters may produce significant scatters in the observed data. 
Clearly, one has to apply the model directly to the observed data to assess its merit. For 
some events with detailed observations over a large dynamical range, one may use the large 
scale turbulence spectrum as input to test the diffusion approximation we have promoted 
here. In some cases, the time dependence of the turbulence spectral evolution may also 
need to be considered, especially in cases with strong anisotropy ob served at large scales 
jMatthaeus et al. Ill990l : basso et al.ll2005l : bsman fc Horburv II2007I ). 



Our current model only consider the Alf ven- Cyclotron branch. To explain the observed 
spectrum above 10 Hz (IDenskat et al. Ill983l ). one may have to include the fast wave mode 
branch, which extends up to the electron cyclotron frequency (Appendix A). Another impor- 
tant observational feature is the high magnetic helicity and the enhancement of the parallel 
magnetic field fluctu a tions in the dissipation ran ge as compared with the inertial range 
(ILeamon et al. 1 1 19981 . Il999l ; lHamilton et al. 1 120081 ) . Considering the polarized characteris- 
tics of the waves, our diffusion model can readily produce results that can be compared with 
these observations. 



7. Conclusion 

The free energy dissipation in collisionless astrophysical plasmas plays crucial roles in 
our understanding of all kinds of nonthermal phenomena. Given the large amount of energy 
inferred from observations of these nonthermal astrophysical sources, the free energy must 
come from large scale structures, while the plasma heating and particle acceleration usually 
start on microscopic scales. An energy cascade from large to small scales is therefore neces- 
sary. Strong turbulence is expected to mediate this energy dissipation in many astrophysical 
sources, and high energy particles, which are responsible for most of the observed emission, 
compete with the low energy background particles to share this energy. There are currently 
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no theories, which can address the energy partition between low and high energy particles 
under different astrophysical conditions. And the major challenges reside in the not-well- 
known nature of turbulence in magnetized collisionless plasmas. Although there is already 
much research both numerical and theoretical on turbulence, none of them gives quanti- 
tative predictions on the nature of the high energy particles produced in this free energy 
dissipation process. In this paper we have proposed the treatment of turbulence cascade 
using the diffusion approximation in the 2D wavevector space, with the aim of treating the 
free energy dissipation mediated by magnetized turbulence self-consistently. To simplify the 
problem, we assume that the low energy particles reach a thermal equilibrium. One there- 
fore can study the turbulence cascade and damping by solving the kinetic equation for the 
turbulence power spectrum numerically. 

The diffusion tensor is the central part of the model. Based on previous analytic and 
numerical results, we have tested two forms for the diffusion tensor; one locally isotropic and 
one anisotropic. We compare the cascade process of Alfven-cyclotron turbulence for these 
two models. As expected the resultant spectra tend to be more anisotropic for the anisotropic 
tensor but in both cases the energy flux is predominately in the direction perpendicular to 
the large scale magnetic field, which is qualitatively (but not quantitatively) similar to the 
Goldreich & Sridhar (1994, 1995) critical balance cascade. The isotropic model for Alfven 
turbulence has a critical balance of kn cx fc 3//4 determined by the di ffusion tensor chosen 



The c ritical balance of the anisotropic tensor is almost identical to the iGoldreich fc Sridhar 



( 119951 ) relation. 



As cascade proceeds to smaller scales (k — > Q p /va), the dispersion relation begins to 
deviate from the simple form valid in the MHD regime. We have included this effect (using 
two reasonable forms for Tyy) and shown that the turbulence spectrum steepens when we 
reach this regime. We also show that the inclusion of thermal damping (based on the hot 
plasma dispersion relation) induces cutoffs in the spectrum as the waves and turbulence 
dissipate and heat the plasma or accelerate particles. These cutoffs are also very anisotropic 
and occur at different wavenumbers for different angles of propagation. As a result of these 
anisotropies, the angle averaged spectra resemble a broken power law both in the wavenumber 
or frequency domain. We emphasize that this broken power law spectrum requires anisotropy 
in the cascade as well as damping. 

We compare these spectra with those observed in the solar wind by Leamon et. al. 
(1999) and find a good fit to the observations with the isotropic diffusion tensor model 
with essentially one free parameter, i.e., the gas density, which is readily observable with the 
current instruments in space. Based on this model, we also make predictions on the variation 
of the break frequency and spectral "index" in the dissipation range with the angle between 
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the magnetic field and solar wind velocity, Alfven Mach number, and plasma beta. Some of 
these results appear to be in agreement with t he weak correlations fou nd recently with the 



statistical analyses of a large sample of events (IMarkovskii et al. 1120081 ) 



Our model is d istinct from existing ID models (IStawicki et al. I l200ll ; iGaltier 1 12006 



Howes et al. 



2008af) for the observed broken power law spectrum of solar wind turbulence. 



Stawicki et al. I (120011 ) essentially used the observed spectrum to infer the re quired d i ffusion 



coefficient for the wave power spectrum. The weak turbulence theory o f IGaltier I (120061 ) 
retains only cascade in the perpendicular direction. In lHowes et al. I (j2008al ) model, a critical 
balance is assumed, which essentially leads to a ID problem. Due to difficulties in producing 
a broken power law spectrum with the dissipation processes, all these models invoke some 
kinetic effects of the wave dispersion on the t urbulence energy cascade. Although for the 



particular event observed by lBale et al. I (120051 ). it appears that kinetic effects play dominant 
roles. The varieties of the observed high frequency solar wind turbulence spectrum suggest 
that dissipation processes should play an important role since the kinetic effect induced 



spectra l break does not vary significantly for different observation intervals. In lHowes et al. 



( J2008ah model, both kinetic and diss ipation processes are considered to p roduce a varieties 
of solar wind turbulence spectrum. ICranmer fc van Ballegooijen I ( 20031 ) were the first to 
study a truly 2D model. However, the model was co nstructed with a much com plicated 
approach to recover the critical balance proposed by iGoldreich fc Sridhar I (119951 ) and it 
hasn't been applied to the solar wind turbulence directly. With the diffusion approximation, 
our model is simplified significantly. Although the diffusion approximation may not be well 
justified, especially for weak turbulence, observations can be used to guide the construct of 
the diffusion tensor and the wave kinetic equations. Direct and more systematic applications 
of these models to the observed solar wind spectra will be able to test them. 
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APPENDICES 



A. Cold Plasma Dispersion Surface 

The stochastic particle acceleration theory is built on wave-particle resonant interac- 
tions, and the plasma turbulence cascade also highly depends on the wave-wave resonances. 
To investigate these resonant interactions, and the acceleration and cascade processes, an 
exact form of the plasma dispersion relation is required. However, a complete treatment of 
the coupled oscillations of particles and fields is too complicated to provide even a stable 
numerical solution that can be integrated into the numerical study of wave-particle or wave- 
wave resonant interactions. By introducing a cold plasma approximation, Stix (1962) solved 
the complete dispersion relation analytically. In this section we discuss the solution of the 
cold plasma dispersion relation and the procedure of isolating individual wave mode for the 
turbulence cascade and particle acceleration study. 

Although the solution of c old plasma disp ersion relation is segmented and has poles at 
particle cyclotron frequencies ( SwansorH ll98ol ) . each physical mode of dispersion surface is 



both continuous and smooth almost everywhere. Thereafter, one can isolate each continuous 
mode with segmented functions that gives k(u>,9) relationship by one-to-one mapping. In this 
subsection, we discuss the dispersion relatio n and the range of each wave mode specifically 



with the notation given by ISwanson I (J1989J). 



Alfven-He Cyclotron branch, 



R and L term reach their first pole (fl a ) and lim k = oo. This cutoff frequency is the end 
of the Alfven branch. 

Fast-Proton Cyclotron branch, 



1 B + F 



B-F 
2A 



Similarly, the fast branch cuts off at Q p with k\\ — > oo. Although Q a is a pole for R 
and L, lim — — — ij m — — ^ 00) [ e the dispersion surface is continuous. The 
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switch of sign at Q a can be easily explained in a simplified cas e where 9 = 0. At 9 = 



F = 2 V P 2 D 2 = 2 | | and the wave solution (ISwanson Ill989l ) simplifies into 

/ I pn I \ 1 / 2 

k = uj I S ± J =lu(S± sign(P D)D) 1/2 (A2) 

Since P < for all a; G [0, 1], and D swiches sign at fi a whereas (P 2 D 2 ) 1 ^ 2 does not, the 

discontinuity is only introduced by the attempt of writing k(u) in an explicit form. As a 

B±F 

result, by switching the sign in — at D = 0, we can follow the continuous dispersion 



surface. 

Whistler branch, 



2A 



k 



B + F r r^^ 



uJ 



B-F 



CO G (Q p , LO e ) 



B + F 

Due to the strong He cycltron (left-handed) damping around fl a , — becomes pure 

imaginary at [fi a , cji), where c^i is the smaller root of equation R(oj) = 0. This id\ is 
the starting point of Whistler branch {u)\ ~ 0.585fi p very insensitive to density and field 
strength). On the other hand, at u e we get P = 0, the formula reach the pole in parallel 
direction (6 = 0), and the Whistler branch cuts off at the electron Langmuir oscillation. 

In the parallel direction, the fast branch starts with electron cyclotron branch and turns 
into proton cyclotron branch; the whistler branch, on the ot her hand, starts with p roton 



cyclotron branch and turns into electron cyclotron branch (IPetrosian fc Liu 1 120041 ) . As 
shown in Figures [T3l and [T4"t this turning point is also the crossing point of electron cyclotron 
and proton cyclotron branches. If we assume turbulence can only start at large scales and 
low frequency, this would be the only point where whistler branch and even higher frequency 
branches obtain energy from the turbulence cascade. This crossing wavenumber k c and the 
efficiency of reverse cascade determine the low wavenumber cutoff on energy spectrum of 
proton cyclotron branch, k m [ n , which in turn determines the acceleration efficiency on ions, 
especially 3 He due to its special charge to mass ratio (see Liu et al. 2006 for a detailed 

discussion of how k m[n affects 3 He spectrum). Thereafter, it is worthwhile to calculate this 

B + F B-F . 
crossing point. By solving — = — (i.e. D = 0) or 

-^l -^l -n a u; 2 a _ 
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we get uj 2 = 0.6150, £3 and the corresponding wavenumber 

k c = ^460.502 - 10- 12 a + 0.377Op/c (A4) 

where a = uj e /^e — y/4:irn e m e c/ B is the only free parameter for the dispersion relation with 
the cold plasma approximation. 




Fig. 13. — The dispersion relation for a cold fully ionized H and He (10% by number) 
plasma at different angles (8 = 0, vr/10, 7r/5, 37r/10, 27r/5, 7r/2 with arrows pointing to the 
direction of increasing angle). In this example, the only free parameter a = uj pfi /Vt e = 0.5. 
Left: Aflven (dashed), Fast (solid) and low frequency segment of Whistler branch (dotted). 
Right: High frequency segment of Whistler branch (dotted), Upper Hybrid (solid) and the 
two EM branches (dashed). 



Upper hybrid branch, 



UJ 



k 



B + F 
2A 



uj G [uj 3 , tt £ 



1 B - F 



2A 

The upper hybrid branch starts at u^, which is the second root of R(u) = 0, and it be- 
comes electron cyclotron wave in parallel direction. It extends to higher frequencies in other 
direction and reaches the resonace frequency u; 4 in perpendicular direction] 15 ! 



14 



This result depends on the particle mass and relative abundance only, which are approximately constant. 

RL 

15 In the perpendicular direction n 2 = so u>4 is the root of S = 
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At even higher frequency, there are two electromagnetic wave (EM) branches. The lower 
one is, 

IB-F 

2A w e [u p>e ,u e ) 



B + F > 



and the higher branch is 



B-F 

k = u\l 2A ^G[u; 5 ,oo), (A5) 



where u)$ is the root of L = 0. 

With these segmented functions, one can construct an analytical one-to-one k(u, 6) 
mapping. Thereafter, the explicit dispersion relation u;(k) can be easily solved numerically. 
Figure [13] shows the calculated dispersion surface for a typical solar flare conditions. The 
accuracy and stability of the dispersion relation obtained by solving these equations will 
facilitate the numerical study on wave particle interaction or turbulence cascading. 



B. Damping and Hot Plasma Dispersion Relation 

The cold plasma approximation ignores any thermal motion of charged particles, and 
therefore, limits the resulting dispersion relation to scales larger than the thermal kinetic gy- 
roradius. In the damping phase, there are strong couplings between the charged background 
particles and turbulent motion, and the background particles are energized. This process is 
a major mechanism for heating the solar flares and other astrophysical plasmas (Petrosian 
et al. 2006, Bittner et. al. 2007). On the other hand, the process also significantly damps 
the plasma turbulence and affects the dispersion relation. To study the dispersion relation 
as well as damping effects (i.e., the imaginary part of wave frequency in the damping 

range, one needs to study the collective motion of thermal particles. It is easier to use the 
distribution function _/j-(r, v, t) to represent the charge and current density, 

P = X> / d 3 vfj (Bl) 

3 J 

J = YtQjf d'vvh (B2) 

3 J 

where the subscript j stands for species of particle. The particle distribution function fj 
satisfies Vlasov equations, which includes, 

~if + v • Vfj + — (E + v x B) • V v f 3 = (B3) 
at rrij 
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and the Maxwell equations. 

There are many linear, quasilinear and nonlinear approximations to solve the Vlasov 
equation set. A detailed introduction to these solutions can be found in Swanson's Plasma 
Waves (1989). To study turbulence cascading and particle acceleration, we only use an an- 
alytical approximation for parallel propagating waves and the numerical nonlinear approxi- 
mation with the WHAMP (Waves in Homogeneous, Anisotropic, Multicomponent Plasmas 
by Ronnmark 1982) code for non-parallel propagating waves. 



The parallel propagating waves are extensively studied in stochastic acceleration theory 
due to its simplicity and efficiency in accelerating particles. The theory is used to explain 
the electron and ion spectra and 3 He enrichment from solar flares (Petrosian & Liu 2004; Liu 
et al. 2004; 2006 on 3 He rich impulsive solar energetic particle events). These studies show 
that the turbulence energy spectrum and its high wavenumber cutoff due to damping are 
important factors that determine the accelerated particle spectra and relative abundance. In 
this subsection, we calculate the damping rate and the cutoff point for parallel propagating 
waves with the linear approximation. 

For parallel and quasi-parallel propagating waves, the Dielectric Tensor can be simplified 
and provide the dispersion relation for right-handed cyclotron wave (i?-wave) and left-handed 
cyclotron wave (L-wave) (Swanson 1982, p. 158), 



B.l. Damping of Parallel Propagating Waves 




uok z Vj 



( 




2 



(B4) 



3 



where u = oo r +iuJi is the complex wave frequency , Vj here is thermal velocity, ( n j = 
is a shorthand, and 




(B5) 



is the Plasma Dispersion Function that can be approximated with 




(B6) 



Except for the electromagnetic waves, Qj 3> kvj for all particle species and all possible 



wavenumb er s 1^1 As a result, the imaginary part of Z(Q for parallel propagating wave is 
small for most frequencies but increase dramatically around particle's gyrofrequency. For 
example, despite the summation of all species of particles in Equation flB4p . the contribution 
from ions to the i?-wave (electron cyclotron wave is the only i?-wave for solar flare or solar 
wind) dispersion relation is far less than that from electrons. Therefore by considering 
electron cyclotron term only, Swanson (1982) obtained the damping rate for the i?-wave 
with the first order approximation (Equation IB6|) . 



U) r 
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pe 



UrkVe 



2 + 



— exp 



UJr 



a 



(B7) 



For the L-wave (for solar flare or solar wind, L-wave includes the helium cyclotron and proton 
cyclotron wave), the contribution from electron cyclotron term can be ignored. However 
due to their close gyrofrequency, the contributions of proton and helium cyclotron term 
become comparable. Thereafter, we generalized Swanson's (1982) derivation to include all 
the particle into the calculation and obtain the damping rate for parallel propagating waves: 



U) r 



E 



uj r kv s 



exp 



UJ r 



kv e 



2/l ^ ^ ^UJy c^i 
UJ^ ' UJ r (uJ r — € S Q S 
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When damping rate overwhelms cascading rate from below, the turbulence energy spec- 
trum cuts off, and the k max is obtained by solving t~J, = = 2ui. Figure [T41 illustrates 
this procedure of finding k max by assuming a Kolmogorov cascading rate (Zhou & Matthaeus 
1990). From the figure We can also see that the cyclotron damping rate grows extremely fast 
(exponential of a square) at particle's cyclotron frequency and thus k max becomes insensitive 
to the cascading rate. Thereafter, by fitting to the numerical results of k max , we am able to 
approximate k max for proton and helium cyclotron branch with plasma temperature and a 
only, 

-0.15 

PC: /,',„„,.: (HOn + 20) ( 1 - 50a + 8 (B9) 



HeC: k max = (50a + 5) 



—) 

10 6 K J 



-0.16 



10a + 6 



(B10) 



16 Although cold plasma approximation suggests that k diverges at particle's gyrofrequency, the damping 
rate at these frequencies diverges too, which prevents the dispersion relation extending to a large wavenumber. 



-46 - 



Dispersion Curve, k mln and k mox (0 = 0) 




10 20 30 40 50 60 

kc/n p 



Fig. 14. — Dispersion relation for parallel propagating waves and its high and low wavenum- 
ber cutoff, k min and k max . Upper panel: The low wavenumber cutoff caused by the topology 
of dispersion curve i.e., the crossing between EC and PC branch. Lower panel: The high 
wavenumber cutoff generated by cyclotron damping overwhelming cascading, we assume a 
simple Kolmogorov cascading rate for both PC and HeC branch in this case. 



Note that although the formula (IB 7ft and (1B8I) gives the damping rate u>i for all k, the 
Lauren series approximation of plasma dispersion function (Equation IB6|) fails at particle's 
gyrofrequencies ( lim ( = 0, where subscript s stands for different species of particles). At 

these ranges u r ~ Q s , one can simplify Equation (1B4[) and show that u>i <C u r (Swanson 
1989). Since the accurate number of uJi becomes unimportant at the range LO{ > u r , we 
simply extrapolate the approximated damping rate with a power law at u r ~ Q s region for 
the turbulence cascade and wave-particle interaction studies. 
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